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by 
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SUMMARY 


During  the  past  30  years  there  has  been  a  resurgence  of  interest  in  the 
classical  orbital  boundary-value  problem  of  Lambert,  largely  because  of  its 
relevance  to  space  rendezvous  and  interception.  The  most  notable  contribution  to 
the  subject  was  by  Lancaster,  Blanchard  and  Devaney,  in  1966,  but  more  recent 
researchers  have  failed  to  build  on  that  work;  the  present  Report  is  aimed  at 
remedying  this  neglect  by  providing  details  of  a  universal  solution  of  Lambert's 
problem  based  on  the  approach  of  Lancaster  et  al.  In  particular,  the  Report 
presents  starting  formulae  for  Halley's  cubic  iteration  process,  used  for  evalua¬ 
tion  of  the  unknown  parameter,  x  ,  at  the  heart  of  the  approach;  this  process 
always  gives  highly  accurate  values  of  x  after  three  iterations. 

A  Fortran-77  computing  procedure  for  a  general  solution  of  Lambert's  pro¬ 
blem  has  been  developed, „aad,.its  three  main  subroutines  are  listed.  Details  are 
given  of  the  testing  oflthis  procedure. 

>  .  . 

Much  of  the  Report  is  devoted  to  a  classification  of  the  set  of  all  Lambert 
problems,  and  to  a  discussion  of  various  geometric  and  physical  aspects. 
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1  INTRODUCTION 

The  'orbital  boundary-value  problem',  constrained  by  two  points  and  an 
elapsed  time,  is  usually  associated  with  the.  name  of  Lambert1,  though  Euler  had 
studied  the  problem  some  20  years  before  Lambert  (but  only  for  parabolic  orbits); 
other  celebrated  mathematicians  whose  names  are  associated  with  the  problem  and 
its  solution  include  Gauss  and  Lagrange,  Thus  it  is  a  problem  of  classical 
relestisl  mechanics,  and  one  that  (like  the  solution  of  Kepler's  equation)  con¬ 
tinues  to  attract  the  attention  of  mathematicians  searching  for  solution  pro¬ 
cedures  of  ever-greatcr  generality,  accuracy  and  efficiency.  Good  text-book 
introductions  u,e  Lo  be  found  in  Kora  2  to  4 ,  whilst  Refs  5  to  26  are  studies, 
chronologically  listed,  from  the  last  30  years;  the  outstanding  paper  on  this 
list,  though  from  as  far  back  as  1966,  i.3  the  one  by  Lancaster,  Blanchard  and 

g 

Devaney  .  Classically,  Lambert's  problem  arose  as  a  core  component  in  the 
determination  of  an  orbit  from  three  observations  of  direction  alone,  the  central 
observation  being  used  (on  a  trial-and-error  basis)  as  a  source  of  the  missing 
distance  data  for  the  other  two  observations.  In  the  Space  Age,  with  direction 
measurement  a  commonplace,  the  solving  of  Lambert's  problem  is  directly  appli¬ 
cable  to  the  important  subject  of  orbital  rendezvous. 

Lambert's  problem  may  be  stated  as  follows:  an  (unperturbed)  orbit,  about 
a  given  inverse-square-law  centre  of  force,  C  say,  is  to  be  found  connecting 
two  given  points,  P^  and  ,  with  a  flight  time  At(“  t^  -  t ^ )  that  has 

been  specified*.  The  problem  muse  always  have  at  least  one  solution  and  the 
actual  number,  which  we  denote  by  N  ,  depends  on  the  geometry  of  the  triangle 
CPjP2  and  the  value  of  At  -  it  is  assumed,  for  convenience  and  with  no  loss 
of  generality,  that  At  >  0  . 

To  get  an  immediate  feel  for  the  problem,  let  us  suppose  first  that  the 
triangle  CP^  is  not  degenerate,  so  that  8  ,  the  angle  subtended  by  P ^ 
at  C  ,  lies  between  0  and  it  .  Then  it  would  appear  there  must  be  at  least  two 
solutions,  since  an  orbital  path  (in  the  plane  CP^  )  can  be  found  that  subtends 
an  angle  2n  -  B  (•/)< 2  going  the  'long  way  aroand')  as  well  as  one  that  subtends 
6  .  We  can  avoid  this  duality,  however,  by  supposing  the  direction  of  motion  to 
be  specified  in  advance,  so  that  the  two  angles  can  be  deemed  to  define  different 
problems.  There  is  a  further  complication,  since  if  At  is  large  enough,  other 
paths  (necessarily  elliptic.-'!)  will  Ije  possible,  each  of  which  includes  a  number 
of  complete  revolutions.  It  turns  out  (and  will  be  apparent  when  Fig  2  is 


*  A  list  of  symbols  is  provided  at  the  end  of  the  Report. 
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introduced)  that  paths  incorporating  a  specific  number  of  complete  revolutions, 
m  say,  normally  occur  in  pairs*  thus  as  At  increases,  N  (for  a  given  triangle 
and  specified  direction  of  motion)  is  an  increasing  odd  integer,  apart  from  being 
even  (instantaneously)  at  each  (critical)  value  of  At  at  which  two  new  solutions 
emerge,  coincident  for  that  precise  value  of  At  .  We  simplify  the  approach  to 
multiple  revolutions  by  extending  the  distinction  between  individual  Lambert 
problems  so  that  'O'  is  regarded  as  an  angle  of  unrestricted  positive  magnitude 
defined  by  the  geometry  of  the  path  and  not  just  by  that  of  the  triangle*  we 
write  0^  for  the  reduced  angle  (such  that  0  <  b  <  2n  )  when  it  is  necessary 
to  discriminate.  (We  do  not  have  to  consider  negative  9  ,  because  to  do  so 
would  imply  negative  At  which  we  have  already  excluded  for  convenience . )  We 
have  effectively  redefined  N  such  that  N  ”  1  if  8  <  2n  ;  if  6  >  2n  ,  on  the 
other  hand,  N  »  0,  1  or  2  ,  depending  an  the  relation  of  At  to  the  appropriate 
critical  value. 

Turning  to  degenerate  triangles,  we  consider  these  on  the  basis  of  the 
unlimited  values  of  0  just  introduced,  so  that  B  is  now  kit  for  some 
integer  k  (>0)  .  Then  if  orbital  paths  exist  that  are  not  rectilinear,  their 
number  must  be  infinite,  since  any  plane  through  the  degenerate  triangle  contains 
valid  paths.  If  we  choose  an  orbital  plane  (as  well  as  the  direction  of  motion) 
arbitrarily,  however,  we  have  N  =  0,  1  or  2  ,  if  k  is  odd,  exactly  as  in  the 
last  paragraph*  this  is  actually  the  simplest  of  all  cases  to  deal  with  in  prac¬ 
tice,  though  the  literature  contains  a  number  of  solution  procedures  that  fail 
here  quite  unnecessarily.  But  there  are  real  difficulties  when  k  is  even 
(«  2m),  associated  with  a  type  of  discontinuity  that  is  described  in  section  3. 

The  effect  of  this  discontinuity  is  that  we  would  like  to  be  able  to  distinguish 
the  angle  (kx)_  ,  which  symbolizes  the  representation  of  6  as  2(m  -  1 ) rr 
plus  a  9  of  2ti  ,  from  (kr)  +  ,  which  symbolizes  its  representation  as  2mn 

plus  a  0  of  zero*  if  this  distinction  (or  an  equivalent  one)  is  not  made,  then 

(for  an  arbitrarily  chosen  orbital  plane)  N  «  1  or  3  if  m  =  1  ,  and 
N  =  0,  1,  2  or  4  if  m  >  1  .  The  orbital  path  has  to  be  rectilinear  (when  k 
is  even)  unless  Pj  and  coincide. 

We  can  now  summarize  the  data  involved  in  the  solution  procedure  to  be 
developed  in  the  present  Report.  The  input  quantities  are  the  constant  u 
(strength  of  the  given  force  centre  at  C  ),  and  r^  (equal  to  CP^  and 

CP,,)  ,  9  (the  unrestricted  angle  pjCP2  ^  ant*  ‘  (We  assurae  u  >  0  ,  but 
there  is  also  a  Lambert  problem  when  n  <  0  ;  if  0  <  x  ,  there  is  then  a  unique 
hyperbolic  solution,  wholly  internal  to  the  triangle.  The  transitional  case, 
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with  jj  *>  0  ,  La  of  course  trivial,  but  even  this  can  be  treated  as  a  Lambert 
problem,)  The  output  quantities,  4N  +  1  in  number,  consist  of  N  itself  and 
N  sets  of  four  quantities,  Pta  V_  (radial  velocity)  and  V_  (transverse 
velocity)  at  both  Pj  and  •  It  is  assumed  here  that  values  of  6  equal 

to  (2tmr)_  and  (2mir)  +  can  be  distinguished,  so  that  N  does  not  exceed  2. 

In  reality,  of  course,  the  real-number  system  does  not  permit  this  distinction, 
though  tliis  is  a  somewhat  academic  point  in  a  computing  procedure  that  can  only 
operate  for  the  finite  set  of  computable  numbers!  more  importantly,  with 
'multipLe  revolutions'  of  a  rectilinear  orbit  the  problem  has  become  completely 
academic  anyway,  since  it  involves  at  least  one  infinite-velocity  'bounce'  off 
the  force  centre.  Nevertheless,  we  shall  find  it  advantageous  to  compute  with 
a  pair  of  quantities,  q  (introduced  in  section  3)  and  m  ,  in  place  of  just 
0  ;  this  avoids  the  academic  difficulties,  and  is  also  more  efficient. 

As  with  Kepler's  equation,  Lambert's  problem  has  no  satisfactory  direct 
solution  -  we  have  to  proceed  by  an  iterative  technique  (trial  and  error)  and 
this  inevitably  dominates  the  solution  procedure  being  developed.  The  following 
issues  then  arise,  and  will  be  discussed  in  successive  sections  of  the  Report! 
first  (in  section  2)  the  choice  of  a  suitable  parameter  of  the  motion  to  use  as 
the  iteration  variable  x  (it  is  sometimes  claimed,  for  example  in  Ref  9,  that 
the  problem  is  inherently  a  two-parameter  problem,  with  simultaneous  iteration 
needed  on  both  parameters,  but  this  claim  is  unwarranted)!  second  (in  section  3) 
the  'direct'  algorithm  that  generates  a  quantity  equivalent  to  At  ,  together 
with  such  of  its  derivatives  as  are  required,  from  x  ,  r1  ,  r,,  and  0  j  third 
(in  section  4)  the  iteration  process,  by  means  of  which  successive  x^  (estimates 
of  x  )  are  computed!  fourth  (in  section  5)  the  starting  formula  (or  formulae) 
for  provision  of  x^  ■,  fifth  (in  section  6)  the  basis  for  the  cessation  of 
iteration  (when  'convergence  is  complete'),  and  the  accuracy  obtained  as  a  result; 
sixth  (in  section  7)  the  formulae  for  computing  the  4N  velocity  components!  and 
lastly  (in  section  8)  the  rationale  behind,  and  results  of,  the  testing  of  the 
solution  procedure. 

2  CHOICE  OF  PARAMETER  FOR  THE  ITERATION  VARIABLE 

2 . 1  Lambert's  theorem  and  the  relations  of  L-congrusnce  and  L-similarity 

Foi  the  iteration  variable,  x  ,  it  is  desirable  to  use  a  quantity  that  is 
a  'Lambert  invariant'  of  the  problem,  if  possible.  To  explain  this  (in  section 
2.2),  we  require  a  preliminary  digression  on  Lambert's  theorem;  as  a  result  of 
this  theorem,  and  the  equivalence  classification  of  triangles  that  it  makes 
possible,  individual  Lambert  problems  can  be  divided  into  equivalence  classes. 
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Lambert's  theorem  is  usually  stated,  with  an  extension  of  the  notation 

already  introduced,  as  follows:  if  c  is  the  chord  (aide  P  P2  )  of  the  triangle 

CP,P2  >  th°n  At  (for  a  connecting  orbital  path)  can  be  expressed  as  a 

(multivalued)  function  of  just  three  quantities  (not  counting  u  ),  vi:. 

r  +  r.  ,  c  and  a  ,  this  last  being  the  semi-major  axis  of  the  path  (taken  us 
i  “  327  2-4 

negative  for  hyperbolic  orbits  *  );  many  text-books  prove  the  theorem  for 

elliptic  orbits,  and  a  general  'minimalist*  proof  has  recently  been  given  by 
28 

Sarnccki  .  Defining  s  as  the  semi-perimeter  of  the  triangle,  so  that 
r |  +  r^  =  2s  -  c  ,  and  noting  that  a  is  equivalent  to  the  total  energy  of  the 
motion  per  unit  mass,  we  can  also  express  At  as  a  function  of  s  ,  c  and 
energy . 

It  follows  from  the  theorem  that  triangles  with  the  same  values  of  s  and 
c  are  equivalent,  from  the  viewpoint  of  the  relation  between  At  and  energy, 
and  the  set  of  all  triangles  CP(P 2  can  thus  be  divided  into  equivalence  classes, 
as  illustrated  in  Fig  1.  Each  class  contains  a  unique  (apart  from  orientation) 
isosceles  triangle  with  r  -  r^  ,  and  the  general  class  (with  0  <  c  <  a, 
illustrated  in  Fig  la)  contains  a  pair  of  degenerate  triangles  such  that  one  of 
the  points  P^  and  P2  lies  between  the  other  point  and  C  ;  a  connecting  orbit 
for  either  of  these  degenerate  triangles  is  necessarily  rectilinear.  Classes 
with  c  =  0  (illustrated  in  Fig  1b)  contain  only  a  single  member  each,  which  is 
simultaneously  isosceles  and  doubly  degenerate.  The  other  extreme  (illustrated 
in  Fig  1e),  occurring  when  c  =  s  ,  is  such  that  the  classes  have  their  widest 
membership,  in  regard  to  the  r2'cj  ratios  possible,  though  all  members  are  now 
degenerate;  each  class  contains  a  pair  of  doubly  degenerate  triangles  such  that 
either  P(  or  coincides  with  C  ,  whilst  the  remaining  (singly  degenerate) 

triangles  (infinite  in  number,  as  in  the  general  case)  all  have  P  and  P^  on 
opposite  sides  of  C  .  Connecting  orbits  for  the  singly  degenerate  triangles  of 
this  extreme  case  cannot  be  rectilinear,  on  the  usual  assumption  that  C  is  a 
point  of  reflexion  (at  infinite  velocity)  for  rectilinear  orbits;  however,  a 
connecting  orbit  for  either  of  the  doubly  degenerate  triangles  is  bound  to  be 
rectilinear  (as  with  singly  degenerate  triangles  in  the  general  case).  For  the 
extreme  classes  with  c  -  s  ,  therefore,  it  is  convenient  to  regard  the  term 
'degenerate'  as  referring  only  to  the  doubly  degenerate  triangles.  If,  further, 
we  cease  to  distinguish  between  the  pair  of  degenerate  triangles  with  <  r^ 
and  r^  >  r2  ,  then  we  can  say,  in  all  cases,  that  an  equivalence  class  contains 
exactly  one  degenerate  triangle. 
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Triangles  in  one  of  the  foregoing  equivalence  classes  may  be  described  as 
'Lambert  congruent',  or  L-congruent  for  brevity,  and  introduction  of  the  concept 
of  congruence  suggests  the  allied  one  of  similarity,  just  as  in  elementary 
geometry.  Thus  two  triangles  may  be  described  as  L*simiiar  whenever  they  have 
the  same  value  of  c/s  ,  this  being  a  dimensionless  quantity.  Though  we  continue, 
in  section  2.2,  to  introduce  Lambert  invariance  on  the  basis  of  L-congruence ,  in 
section  3  we  shall  find  that  thinking  in  terms  of  L-similarity,  with  its  wider 
equivalence  classes,  has  the  effect  of  reducing  (by  one)  the  number  of  arguments 
in  the  algorithm  for  the  flight  time. 

2 . 2  Lambert-invariant  parameters,  Lambert-equivalent  problems,  and  the 

parameter  x 

Since  At  is  a  function  only  of  s  ,  c  and  energy,  it  follows  that  the 
equivalence  of  L-congruent  triangles  provides  the  basis  for  a  classification  of 
individual  Lambert  problems  into  their  own  equivalence  classes,  each  such  class 
being  defined  by  the  underlying  class  of  triangles  and  the  given  value  of  At  . 
Then  a  Lambert-invariant  parameter  may  be  defined  as  one  that  has  the  same  value 
for  all  members  of  an  equivalence  class  of  problems.  It  is  unfortunate  that., 
though  s  and  c  are  Lambert<*invariant ,  6  is  not,  which  at  first  sight 
negates  the  virtues  of  the  unrestricted  angle  that  were  noted  in  section  1 .  We 
can  get  the  best  of  both  worlds,  however,  by  taking  0  (instead  of  B)  as  a 
parameter  of  the  general  triangle,  where  0  is  defined  as  being  the  6  for  the 
equivalent  isosceles  triangle:  then  0  can  be  regarded  (like  6)  as  an  angle 

of  unrestricted  magnitude,  (The  quantity  J©  ,  denoted  by  f  ,  was  recognized 

1 9 

as  an  important  parameter  in  the  paper  by  Batttn,  Fill  and  Shepperd.) 

The  energy-equivalent  orbital  parameter  a  is  certainly  Lambert-invariant, 
hut  this  is  not  true  of  e  (eccentricity)  or  p  (semi-latus  rectum) .  The  use 
of  p  as  iteration  variable  is  intuitively  appealing,  because  of  its  direct 
relation  to  true  anomaly  and  hence  to  9  ,  and  it  is  recommended  in  a  paper  as 
recent  as  Ref  24.  But  p  fails,  in  a  somewhat  paradoxical  fashion,  when  9  =  it  . 
The  paradox  is  that  p  is  given,  as  +  r,,)  ,  without  the  need  to 

iterate  at  all  in  these  circumstances:  because  At  is  not  involved  in  this 
formula,  however,  no  further  progress  can  then  be  made  without  iterating  on 
some  other  variable. 

The  advantage  in  using  a  Lambert-invariant  parameter  as  the  iteration 
variable  is  that  its  determination  is  a  numerically  identical  procedure  for  all 
the  individual  problems  of  an  equivalence  class.  The  resulting  'reduction  in 
cases'  is  a  very  practical  consideration  for  the  solution  procedure  to  be 
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developed,  not  least  when  it  comes  to  the  number  of  tests  that  have  to  be  run. 
Thus  a  Is  an  immediate  candidate  for  iteration  variable,  but  its  direct  use 
would  be  unsatisfactory  because  (as  shown  in  Ref  2,  for  example)  orbital  paths 
with  a  particular  value  of  a  occur  in  pairs  or  not  at  all.  There  is,  in  fact, 
an  upper  limit  (corresponding  to  minimum  possible  energy)  to  the  possible  values 
of  1 /a  ;  it  is  given  by  2/s  ,  for  which  value  the  pair  of  paths  coincide,  it 
follows  from  this  that  if  we  write 


then  x  is  a  satisfactory  substitute  for  a  ,  such  that  the  choice  of  sign  in 
the  implied  square  root  distinguishes  between  the  two  paths  of  each  pair.  The 
parameter  x  is  universal  (defined  independently  of  the  type  of  orbital  path)  , 

unlike,  for  example,  the  parameter  used  by  Sun,  Vinh  and  Chern  in  their  recent 
26 

paper  ;  it  also  has  the  advantage  of  being  non-dimensional,  which  facilitates 
the  switch  from  L-congruence  to  L-simiiarity . 

The  parameter  x  ,  as  just  introduced,  is  the  iteration  variable  used  in 

Q 

che  milestone  paper  by  Lancaster,  Blanchard  and  Devancy  (of  which  Ref  10  is  a 

somewhat  expanded  version)  underlying  much  of  the  work  reported  here.  It  has 
.  28 

been  shown  by  Sarnecki  chat  x  has  a  dynamical  interpretation,  being  a  non- 
diraensionalized  value  of  the  velocity  in  the  (rectilinear)  solution  of  the 
Lambert  problem  for  the  degenerate  triangle  that  is  L-congrucnt  to  the  given  one 
and  is  such  that  Tj  i  >  x  is  positive  or  negative  according  to  whether  the 
direction  of  V  (now  pure  radial)  is  inward  or  outward,  and  (in  principle  -  see 
also  section  3)  this  resolves  the  ambiguity  in  the  sign  of  x  outstanding  from 
the  last  paragraph.  (Battin  finds  an  interpretation  for  x  in  terms  of  the 
actual  problem,  whilst  a  geometrical  interpretation  of  the  parameters  of  the 
classical  Lambert-Euler  equation,  also  involving  the  rectil inearly  equivalent 
problem,  is  given  in  Ref  20.)  Sarnecki's  interpretation,  which  requires  P^  to 
be  the  more  distant  point,  involves  an  intrinsic  (and  perhaps  surprising)  lack 
of  symmetry;  thus  s  =  max  fr^,  r^)  when  Pj  and  are  connected  by  a 

rectilinear  orbit.  It  is  clear  from  equation  (1),  now  that  there  is  no  longer 
an  ambiguity  in  the  sign  of  x  ,  that  |x|  <  1  for  elliptic  orbits,  x  -  1  for 
parabolas  and  x  >  1  for  hyperbolas.  Values  of  x  <  -I  do  not  arise;  they 
would  be  associated  with  negative  values  of  at  .  There  is  an  apparent  distinc¬ 
tion  (asymmetry)  between  the  elliptic  and  hyperbolic  paths  corresponding  to  a 
given  pair  of  points  P^  and  Pj  :  if  a  positive  value  of  1/a  ,  less  than  2/s  , 
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is  given,  equation  (1)  gives  two  values  of  x  ,  numerically  equal  but  opposite 
in  sign,  and  to  each  value  there  corresponds  an  elliptic  orbital  path;  if  I /a 
is  negative,  on  the  other  hand,  only  the  positive  root  of  (1)  is  legitimate, 
implying  only  one  hyperbolic  path.  But  suppose  that,  for  a  given  value  of  6 
less  than  u  ,  we  consider  the  path  for  2tr  -  0  as  well  as  the  path  for  0  ; 
then  for  negative  1 /a  there  are  two  hyperbolic  paths  that  are  quite  distinct 
(though  having  the  same  value  of  x) ,  whereas  tor  positive  1 /a  the  two 
'additional'  elliptic  paths,  one  associated  with  positive  x  and  one  with 
negative  x  ,  are  merely  the  'orbital  complements'  of  the  two  original  paths 
(with  reversed  signs  of  x) .  Thus  the  apparent  asymmetry  between  the  two  types 
of  orbit  (and  parabolic  orbits  behave  in  the  same  way  as  hyperbolic  orbits)  has 
no  real  significance.  The  nature  of  the  various  orbital  paths  is  well  illus¬ 
trated  by  Figs  3.7  and  3.8  in  chapter  3  of  Ref  2. 

3  THE  ALGORITHM  FOR  COMPUTING  At 

An  algorithm  for  At  ,  in  terms  of  x  ,  requires  a  pair  of  Lambert- 

invariant  parameters  to  specify  the  relevant  triangle.  By  thinking  in  terms  of 

L-aimilarity ,  rather  than  L-congruence,  however,  we  can  reduce  the  number  of 

triangle  parameters  from  two  to  one,  so  long  as  the  output  is  made  non- 

8  1 0 

dimensional;  we  follow  Lancaster  et  al  ’  in  replacing  t  by  T  ,  where 


The  choice  of  s  ,  as  'length  scale'  in  equation  (2),  pays  dividends  (and  was 

also  made  in  some  of  the  papers  by  Battin  and  his  colleagues,  ey  in  Ref  18).  In 

the  classical  procedure  of  Gauss,  on  the  other  hand,  the  length  scale  is 

effectively  Fs(s  -  c7  ,  whilst  in  the  modification  of  Gauss's  method  due  to 
23 

Battin  and  Vaughan  it  is  somewhat  more  complicated  -  Ref  23  is  perhaps  the  most 
interesting  of  the  papers  since  the  pair  by  Lancaster  et  al,  and  Appendix  A  is 
devoted  to  a  brief  discussion  of  this  Report  and  the  underlying  method  of  Gauss. 
Finally,  in  Ref  26  the  scale  length  is  effectively  2s  -  c  ,  ie  +  r^  .  The 
advantage  of  s  over  any  of  these  alternatives  is  that  T  (as  opposed  to  the 
quantities  corresponding  to  it)  is  monotonic  with  respect  to  0  (for  fixed  x) , 
which  makes  for  the  essential  simplicity  of  Fig  2,  introduced  in  the  next  para¬ 
graph;  further,  we  have  already  seen  that  s  is  the  quantity  that  is  directly 
associated  with  the  minimum-energy  path.  The  importance  of  using  s  as  the 
Length  scale  must  not  be  exaggerated,  however,  as  the  use  of  one  of  the  alterna¬ 
tives  does  not  affect  the  progress  of  the  iteration  process,  except  of  course  if 
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the  male  becomes  zero,  in  particular  with  Gauss's  scale  when  s  =  c  .  (The 
presence  of  p  in  (2)  debars  the  solution  of  Lambert's  problem  when  p  is  zero; 
x  is  then  infinite,  since  a  “  0  for  the  hyperbola  into  which  the  solution 
degenerates,  and  we  can  include  this  case  in  a  general  algorithm  if  we  sacrifice 
full  Lambcrt-Lnvariance  and  solve  for  xfjT  rather  than  x  .) 

The  quantity  T  is  a  function  of  the  two  parameters  x  and  0  (defined 
in  section  2.2),  but  Lancaster  et  aL  use  q  ,  rather  than  O  ,  where 

q  =  cos  !  9  /(1  +  sin  i  0^)  =  tan  i  (n  -  Qr)  .  (3) 

here  0^  is  l)  reduced  to  the  range  (0,  2n  1 ,  so  that 

sini0r  “  *  (4) 

Clearly,  (3)  and  (4)  lead  to  the  simple  result  that 

=  ( 1  -  sin  J  0^) /( 1  +  sin  i  0  )  =  1  “  c/s  ,  (3) 

but  (5)  does  not  define  the  sign  of  q  ,  whereas  (3)  does.  We  also  have,  from 

(3)  and  (4), 


q  =  (1-c/2s)  cos  J  0^ 


(6a) 


which  is  just  a  particular  case  of  the  result  that,  for  the  general  triangle, 
follows  from  the  standard  'cosine  formuLa1  and  may  be  expressed  as 


q 


cos  i  0^ 


(6b) 


we  refer  to  (6a)  and  (6b)  in  Appendix  A. 

When  q  is  the  second  parameter  of  the  T-function,  m  is  required  in 

addition  (as  a  third  parameter)  to  specify  the  number  of  completed  circuits. 

This  brings  us  to  Fig  2,  which  plots  (as  in  Refs  8  and  10)  T  against  x  for 

particular  values  of  q  and  m  (corresponding  to  selected  values  of  0) .  The 

relation  (almost  linear)  between  q  and  0^  is  indicated  by  Fig  3,  which  also 

plots  0^  Q  and  0^  m  against  0r  .  Here  0^  Q  is  a  quantity  that  may  be 

regarded  as  a  'time-linearized'  version  of  0  ,  ie  it  is  defined  over  the  range 

(0  ,  2 u)  in  such  a  way  as  to  make  T  a  linear  function  of  G_  _  when  x  =  0  ; 

1 
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similarly,  0^,  rn  becomes  time-linearised  as  x  tends  to  »  .  These  two  curves 
are  much  less  linear,  the  non-linearity  for  0^  g  being  pronounced.  (As 
x  ■*  ”  ,  T  behaves  like  2x_,(1  -  q|q|)  ,  so  0T  ^  =  v()  -  q|q|)  .) 

The  most  striking  feature  of  Fig  2  consists  in  the  gaps  (unrealizable 
regions)  that  occur  in  the  part  of  the  figure  associated  with  elliptic  orbits. 
Thus,  as  0  increases  through  a  multiple  of  2n  ,  for  a  fixed  value  of  x  3uch 
that  0  <  j x j  <  1  ,  q  jumps  from  -1  back  to  +1  and  there  is  a  jump  in  the  value 
of  T  .  This  brings  us  to  the  difficulty,  referred  to  in  section  1  in  relation 
to  the  use  of  0  ,  in  the  functional  representation  of  T  as  T(x,0)  rather 
than  T(x,q,m)  j  when  0  »  2mi7  and  x  ^  0  ,  in  fact,  T(x,0)  is  not  unique, 
the  values  of  T(x,  -1,  m  -  1)  and  T(x,l,m)  not  being  the  same.  Though  this 
is  not  a  major  defect  in  the  T(x,0)  representation,  as  we  saw,  it  combines  with 
the  greater  efficiency  of  the  T(x,q,m)  representation  to  make  us  regard  the 
latter  as  'standard1. 

We  can  account  for  the  unrealizable  regions  of  Fig  2  by  studying  (for  a 
fixed  x  with  |x|  <  1  )  the  variation  of  T  with  0  for  the  isosceles 
representatives  of  the  equivalence  classes  of  triangles.  If  x  =  0  ,  the  con¬ 
necting  orbital  paths  are  all  minimum-energy  trajectories,  by  equation  (1),  and 
there  is  no  jump  at  9  “  2mir  (we  now  write  9  ,  rather  than  0  ,  since  they  are 
identical  for  isosceles  triangles) .  When  x  ^  0  ,  for  which  a  number  of  cases 
(with  a  fixed  value  of  jxj  ,  vis  0.5)  are  illustrated  in  Fig  A,  consider  first 
the  situation  for  9=0,  regarded  as  the  limit  case  as  6  reduces  to  zero. 

Then  our  triangle  is  degenerate  as  well  as  isosceles,  but  the  orbital  path  does 
not  itself  have  to  be  degenerate  (rectilinear).  In  fact,  any  orbit  through  P^ 
trivially  passes  through  P^  after  zero  time  (since  the  points  coincide),  and 
this  is  reflected  in  the  zero  value  of  T  when  x  >  0  (left-hand  illustration 
of  Fig  4a);  uniqueness  here  derives  from  the  isosceles-triangle  assumption,  the 
limit  case  of  this  being  such  that  che  velocity  is  pure  transverse  and  hence  the 

orbit  non-degenerate  -  for  x  =  1//2  ,  indeed,  the  orbit  is  circular.  (The 

2 

general  formula  for  eccentricity,  under  these  conditions,  is  |l  -  2x  |  .)  An 
entirely  different  path  is  available,  however,  that  can  be  seen  to  be  intrinsi¬ 
cally  unique.  It  is  given  by  a  degenerate  (rectilinear)  orbit  with  V  directed 
outward  from  ,  such  that  (still  supposed  coincident  with  P^)  is 

reached  after  non-zero  T  ,  even  though  9  -  0  ;  this  applies  when  x  <  0 
(right-hand  illustration  of  Fig  4b),  the  velocity  now  being  pure  radial. 

Now  consider  what  happens  as  9  increases  from  zero  to  2tt  ,  with  a  fixed 
x  that  is  either  positive  or  negative.  The  two  orbital  paths  are  always  such 
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that  T(x  <  0)  exceeds  'C(x  >  0)  ;  and  when  6  s  ti  (Fig  4c) ,  the  path  for 
T(x  <  0)  is  the  reflexion  in  the  line  of  the  orbital  complement  of  the 

path  for  T(x  >  0)  ,  such  that  the  point  half  way  along  the  path  is  perioentre 
if  x  >  0  but  apoeentre  if  x  <  0  .  For  0  >  r  ,  it  becomes  clear,  from 
symmetry  considerations,  that  the  pairs  of  orbital  paths  are  always  shaped  as 
the  orbital  complements  of  the  paths  for  0  <  u  ,  but  with  the  sign  of  x 
reversed  in  the  correspondence  (compare  the  pair  of  illustrations  of  Fig  4d  with 
the  pair  of  Fig  4b).  As  0  comes  up  to  2 ti  (Fig  4u)  ,  the  position  is  as 
follows i  for  x  <  0  ,  the  orbital  path  approaches  the  complete  circuit  of  the 
(non-degenerate)  orbit  from  which  originated  the  infinitesimal  arc  that  applied 
to  0=0  and  x  >  0  ;  for  x  >  0  ,  on  the  other  hand,  the  path  approaches  the 
incomplete  rectilinear  orbit  Lhat  complements  the  path  applying  to  0=0  and 
x  <  0  .  In  the  latter  case,  the  (initial)  radial  velocity  now  has  to  be  inwaras 
and  (in  the  limit)  the  starting  point  is  reached  after  reflexion  (infinite-V 
’bounce1)  from  the  force  centre. 

In  studying  T(x,Q)  beyond  0  =  2n  ,  still  assuming  isosceles  triangles, 
we  have  a  choice  in  the  way  we  ascribe  the  sign  of  x  .  To  follow  Lancaster  at  at 
(and  also  Sarneckl  ),  we  effectively  start  again  with  a  complete  revolution 

behind  us,  so  that  the  position  represented  by  Fig  4e  is  followed  by  the  position 

represented  by  Fig  4a  with  an  orbit  grafted  on.  This  implies  a  double  disconti¬ 
nuity,  so  that  for  both  x  >  0  and  x  <  0  there  is  a  jump  in  the  value  of  T 
at  0  =  2tt  ,  by  AT  say,  equal  to  the  value  that  T  has  when  x  <  0  and  0=0 
(and  hence  equal  to  twice  the  journey  time  from  P^  to  apoeentre).  With  the 
double  discontinuity  associated  with  this  specification  of  the  sign  of  x  ,  we 
have  accounted  for  the  unrealizable  regions  of  Fig  2,  as  promised;  at  0  =  2n 
itself,  there  are  three  (not  four)  possible  values  of  T  ( of  the  original  dis¬ 
course  on  Lambert  solutions  in  section  1)  corresponding  to  a  given  |xj  ,  since 
the  limiting  post-jump  value  for  x  >  0  is  equal  to  the  limiting  pre-jump  value 
for  x  <  0  (same  non-degenerate  orbital  path) ,  both  being  equal  to  the  orbital 
period.  But  this  very  fact  leads  to  the  alternative  way  of  defining  x  when  O 

lies  between  2mn  and  2  (m  +  1)it  for  any  odd  value  of  m  ;  this  simply  reverses 

the  sign  assumed  by  Lancaster  el  al .  Thus,  instead  of  just  ’starting  again’  at 
9  =  2n  ,  we  ’switch  over1  (between  left-hand  and  right-hand  illustrations)  in 
returning  from  Fig  4e  to  Fig  4a.  So  there  is  now  no  jump  in  T  when  x  <  0  ; 
for  x  >  0  ,  on  the  other  hand,  there  is  a  jump  of  2AT  .  When  we  get  to 
0  =  4ir  ,  the  position  is  the  exact  opposite:  there  is  a  jump  for  x  <  0  .  with 
x  as  redefined,  but  no  jump  for  x  >  0  ;  the  sign  rf  x  now  reverts  to  being 
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the  same  as  in  the  definition  by  Lancaster  ct  al.  This  alternation  continues, 
ad  infinitum ,  for  successive  values  of  m  . 

Thus,  if  the  alternative  definition  of  (the  sign  of)  x  is  adopted,  it  is 

possible  to  plot  a  figure  with  only  half  the  number  of  unrealizable  regions,  and 

this  has  been  done  in  Fig  5j  when  the  sign  of  x  is  reversed,  the  sign  of  q 

has  to  be  reversed  as  well,  which  makes  q  always  continuous  with  respect  of  0 

(though  dq/d0  will  have  discontinuities  -  see  again  Fig  3) .  It  is  emphasized 

that  both  methods  of  definition  of  x  are  legitimate;  tne  original  definition, 

due  to  Lancaster  at  al,  leads  to  the  simpler  algorithm,  however,  and  will  be 

adopted  for  the  rest  of  this  Report.  In  this  context,  it  is  worth  remarking 
28 

that  Sarnecki"  effectively  presents  an  understanding  of  T(x,0)  ,  obtained  here 
via  the  isosceles-triangle  representatives  of  the  equivalence  classes,  by  con¬ 
sidering  the  degenerate-triangle  representatives,  Then  increasing  values  of  0  , 
for  a  fixed  value  of  x  ,  can  be  interpreted  in  terms  of  the  'time-line'  of  a 
single  orbital  path,  so  long  as  (now  evolving)  remains  less  than  r^  (fixed 

initial  value),  so  that  x  (and  hence  s)  remains  constant.  From  Sarnecki 's 
viewpoint,  the  unrealized  values  of  T  correspond  to  a  breakdown  of  the  con¬ 
straint  (that  r^  i  r^) ,  but,  regardless  of  this  breakdown,  the  original  defi¬ 
nition  of  x.  cleariy  fits  this  interpretation  more  naturally  than  the  alternative 
one.  In  Ref  23  (the  paper  that  is  discussed  in  Appendix  A),  on  the  other  hand, 
Battin  and  Vaughan  make  explicit  use  of  the  isosceles-triangle  equivalent  of  a 
given  problem,  in  the  solution  that  they  offer  as  an  improvement  on  Gauss's 
method.  (it  is  noted,  for  completeness,  that  Figs  2  and  5  can  be  extended  to 
cover  negative  values  of  T  ,  by  reflecting  the  existing  contours  in  the  origin; 
for  T  <  0  it  would  be  values  of  x  i-  1  ,  not  x  <  -1  ,  that  would  be  impossible, 
however,  and  this  complication  is  one  reason  for  restricting  consideration  to 
positive  At  .) 

Further  light  is  thrown  on  the  evolution  of  the  orbital  paths  for  isosceles 

triangles,  if  the  variation  of  the  eccentricity  is  considered,  and  Fig  6  gives 

contours  of  c  in  the  (x,0)-plane,  with  x  as  'alternatively  defined'.  (For 

2 

0=0  and  x  >  0  ,  we  have  e  =  | 1  -  2x  |  ,  as  already  noted.) 

The  basic  formulae  for  computing  T  from  x,  q  and  ra  were  given  in 
the  papers  by  Lancaster  et  al  and  are  repeated  here  (without  proof)  in  Appendix  B, 
which  gives  some  details  of  the  Fortran-77  subroutine,  TLAMB,  that  has  been 
written  as  the  core  of  a  universal  procedure  for  solving  Lambert's  problem; 
further  details  of  TLAMB  may  be  elicited  from  the  listing  in  Appendix  C.  Fcr 
tlie  present,  a  few  remarks  will  suffice.  The  first  is  that,  in  addition  to  T  , 
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the  subroutine  can  generate,  if  needed,  (partial)  derivatives  (up  to  the  third) 
with  respect  to  x  -  we  denote  these  derivatives  by  T*  ,  T"  and  T1"  . 
Secondly,  the  basic  formulae  of  Lancaster  ct  al  do  not  give  full  accuracy  in  all 
circumstances;  it  should  be  clear  from  Appendices  B  and  C,  however,  that  rounding 
error  is  minimized  in  TLAMB.  Thirdly,  one  of  the  quantities  used  by  the  sub¬ 
routine  is  conveniently  introduced  here,  as  it  is  needed  in  section  7.  The 
quantity  is  z  ,  defined  by 

z  =  +  >4 1  -  q2 * 4  +  q2x2)  .  (7) 

It  has  a  dynamical  interpretation,  similar  to  that  given  for  x  in  section  2, 
following  Ref  28;  thus  z/q  is  a  non-dimensionalized  vaLue  of  the  velocity  at 
the  dollar  point  in  the  degenerate  Lambert  problem  (z  is  never  negative  because 
the  direction  of  this  velocity  is  looked  after  by  the  sign  of  q  ).  The  com¬ 
putation  of  z  by  equation  (7)  is  itself  an  example  of  the  potential  loss  of 

2 

accuracy,  since  in  many  situations  it  is  desirable  that  1  -  q  be  regarded  as 

a  quantity  available  independently  of  q  ;  this  is  true,  in  particular,  when 

in  =  0  and  q  is  close  to  1,  ie  when  (.-)  is  small.  Because  of  this  danger, 

2 

TLAMB  has  1  ••  q  as  an  extra  argument,  the  assumption  being  that  it  may  be 
computed  directly  from  (•)  via  the  formula 

2 

1  -  q  =  c/s  =  2  sin  }  0  /(1  +  sinJO^)  ,  (8) 

which  is  immediate  from  (5). 

4  ITERATION  PROCESS 

Though  the  starting  value,  x^  ,  has  to  be  available  before  numerical 
iteration  can  commence,  the  iteration  process  is  considered  now,  before  we  look  • 
at  starting  formulae,  because  the  iteration  process  drives  the  starting  formulae 
rather  than  vice  verua.  Thus,  if  just  the  basic  Newton-Raphson  method  of 
iteration  is  used,  the  devising  of  starting  formulae  to  cover  all  cases  becomes 

an  almost  impossible  task.  When  the  Halley  process,  found  to  work  extremely 

.  29  30 

well  in  the  solution  of  Kepler's  equation  (for  hyperbolic  as  well  as  elliptic 

orbits),  is  used  instead,  however,  the  task  is  greatly  eased.  This  is  well 

illustrated  by  the  solution  for  x  when  its  true  value  is  0.5  and  0  is  small, 

say  10  Then  (with  the  starting  formula  given  by  equation  (10)  in  the  next 

section)  Xq  is  roughly  double  the  true  x  ,  after  which  the  Newton-Raphson 

process  leads  to  an  X|  very  close  to  zero  and  the  iterative  process  effectively 
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stagnates;  the  Halley  process,  on  the  other  hand,  gives  an  x^  very  close  to 
the  true  value,  after  which  convergence  is  rapid.  (The  reason  for  the  accurate 
Xj  computation,  in  circumstances  such  as  these,  is  that  for  very  Btnall  0  and 
x  »  v'B  ,  T  behaves  like  1/x  ,  which  is  a  bilinear  function,  and  the  Halley 
method  gives  a.,  immediately-correet  solution  for  bilinear  functions  .)  In 
defence  of  the  Newton-Raphson  process,  on  the  other  hand,  it  i9  remarked  that 
it  works  very  well  when  the  value  of  0  approaches  2r  and  |x|  is  small,  so 
long  as  the  appropriate  starting  formula  from  section  5.2  is  used;  this  is 
notwithstanding  the  explicit  warning  (against  using  the  process  in  these  circum¬ 
stances)  issued  at  the  end  of  Ref  10,  an  unwarranted  warning  to  which  there  will 
be  further  reference  in  the  present  paper. 

The  Halley  method  is  essentially  the  Newton-Raphson  method  extended  to  give 
third-order  convergence,  ao  it  requires  T"  (the  second  derivative  of  T  with 
respect  to  x  )  as  well  as  T'  (the  first).  Since  Halley's  method  was  also 
adopted  for  the  iteration  involved  in  a  subsidiary  problem  that  arises  when 
m  s  0  (see  section  5.3),  to  satisfy  an  equation  expressed  in  terms  of  X'  ,  we 
sometimes  also  need  T1"  ;  this  is  why  the  subroutine  TLAMB  generates  derivatives 
up  to  the  third. 

In  the  present  study,  the  iteration  process  is  incorporated  in  the 

Forcran-77  subroutine,  XLAMB,  that  generates  solutions  for  x  and  is  listed  in 

2 

Appendix  D.  The  input  of  XLAMB  consists  of  m  ,  q  ,  1  -  q  (supplied  separ¬ 
ately  for  the  same  reason  as  in  TLAMB)  and  T  ,  and  its  output  is  as  follows: 
the  integer  N  (defined  in  section  1)  that  specifies  the  actual  number  of  solu¬ 
tions  (this  should  get  set  to  0,  1  or  2,  though  a  value  of  -1  is  also  theoreti¬ 
cally  possible,  constituting  the  flag  to  be  defined  in  section  6);  x  ,  a  solution 
when  at  lease  one  exists;  and  x  ,  the  second  solution  when  there  are  two  (x 
is  actually  the  first  of  the  two  solutions  to  be  described  in  section  5.3,  so  it 
is  always  positive) . 

5  STARTING  FORMULAE 
5 . 1  Introductory  remarks 

The  requirement  in  regard  to  starting  formulae  fer  the  iteration  process 
is  for  an  approximation  to  the  inverse  of  the  function  that  (for  given  m  and 
q  )  generates  T  from  x  .  When  m  “  0  ,  there  is  formally  no  difficulty,  since 
there  is  a  unique  x  to  which  the  'starter'  is  to  approximate  (in  view  of  the 
assumption  T  t  0  ).  When  m  >  0  ,  on  the  other  hand,  this  uniqueness  does  not 
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normally  apply:  the  normal  oases  are  of  'no  solution'  and  'two  solutions',  but 
the  case  of  just  one  solution,  arising  as  two  solutions  merge,  is  also  covered. 

All  the  starters  have  been  incorporated  in  the  subroutine  XLAMB. 

5.2  Single-revolution  starters 

For  m  s  0  ,  the  curves  in  Fig  2  are  strictly  monotonic,  except  when  q  ■  1 
and  x  >  0  (excluded  case  of  T  s  0  ) .  As  an  obvious  first  move,  we  can  deter¬ 
mine  the  sign  of  x  by  evaluating  (by  TLAMB)  the  value  of  T  ,  T  say*, 

corresponding  to  x  =  0  .  Then  the  sign  of  x  is  the  sign  of  T  -  T  ,  and  it 

is  used  to  distinguish  two  cases. 

Suppose  first  that  x  >  0  .  Then  we  can  approximate  to  the  contour 
T(x,q,0)  ,  for  the  given  q  (regarded  as  fixed),  by  the  bilinear  curve 

T  =  T2/(T  +  4x)  .  (9) 

O  0 

The  rationale  for  this,  as  a  particularization  of  the  general  bilinear  expression 
(a  +  bx)/(c  +  dx)  ,  is  that  it  satisfies  the  following  three  constraints:  first 
that  T  tends  to  0  as  x  tends  to  infinity;  secondly,  that  T  ■  Tq  for  x  »  0  ; 
and  finally,  that  T'  «  -4  for  x  -  0  .  (T^  has  a  fixed  value,  though  this 

is  not  obvious  from  Fig  2;  when  | q  j  =  1  ,  it  is  not  defined,  but  the  one-sided 

derivatives  both  exist,  being  -8  on  one  side  and  0  on  the  other  so  that  Btill 

has  a  conventional  value  of  -4.)  The  merit  of  the  bilinear  approximation  is  that 
it  is  immediately  invertible,  the  inverse  function  being  also  bilinear;  thus  from 
(9)  we  get,  as  our  starter  for  x  >  0  , 

x  .  T  (T  -  T)/4T  .  (10) 

U  o  o 


If  x  <  0  ,  we  proceed  in  the  same  way,  this  time  approximating  the 
T(x,q,0)  contour  by  the  bilinear  curve  that  is  most  naturally  expressed  as 


T 


T 

o 


4x 

x  +  1 


(11) 


The  first  of  the  three  constraints  that  lead  to  (11)  is  chat  T  tends  to 
infinity  as  x  tends  to  -1;  the  other  two  constraints  are  the  same  as  for  the 
case  x  >  0  .  On  inverting,  to  get  x  as  a  bilinear  function  of  T  ,  we  have 


*  The  zero  suffix  in  T0  is  notated  differently  from  the  zero  suffix  in  x^  ,  to 
reduce  possible  confusion  between  the  two  different  meanings.  1 


02/ 


1? 


T  -  T 

x  .  .  — _ S _ 

01  T  -  T  +  4 

o 


(12) 


This  formula,  (12),  is  as  elementary  as  the  complementary  one,  (10),  but  it  was 
found  to  work  much  leas  well  and  it  was  necessary  to  patch  it;  this  is  why  the 
left-hand  side  of  (12)  has  been  written  x^  rather  than  xQ  .  The  patching  was 
of  a  complicated  nature,  with  some  arbitrary  features,  but  it  will  be  summarized 
(in  the  next  three  paragraphs)  for  completeness. 


Equation  (12)  actually  requires  two  patches.  The  necessity  for  the  first 
became  apparent  for  values  of  0  greater  than  about  1.999v,  and  it  was  ascribed 
to  the  proximity  of  the  (left-hand  half  of  the)  curve  for  0  =  2u  (q  =  -1),  for 
which  the  value  of  T^  should  actually  be  zero  (as  opposed  to  -4,  for 
0  <  0  <  2ir) .  Now  is  automatically  zero  for  a  curve  that  i3  bilinear  in 

x2  (rather  than  x  ) .  Hence  we  are  led  to  consider  the  alternative 
approximation 


T  -  T  (1  +  ix2)/(1  -  x2) 
o 


(13) 
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in  which  only  the  term  Jx  requires  explanation  -  it  gives  T^  its  correct 

value  (bit)  for  0  ■  2it  ,  where  T  =*  2x  .  The  behaviour  of  T"  follows  from 

o 

equation  (B— 16)  of  Appendix  B.  Thus,  if  |q]  <  1  ,  T"  is  defined  and  finite 
for  all  x  ,  but  it  is  unbounded  as  [q|  approaches  unity  and  x  approaches 
zero.  If  |q|  is  actually  constrained  to  be  unity,  however,  T"  is  bounded  in 
the  neighbourhood  of  x  =  0  !  it  is  strictly  undefined  at  T  =  0  ,  but  may  be 
regarded  as  'effectively  defined'  since  it  has  the  same  limiting  value,  equal 
to  3Tq  ,  for  an  approach  to  the  limit  from  either  side.  The  inversion  of  (13) 
may  be  written 

*02  =  “  ^(T  “  V/(T  +  »Vl  '  (1A) 


At  this  stage  we  have  two  possible  formulae  for  our  desired  starter,  the 
first  (and  simpler)  being  normally  to  be  preferred!  thus  the  possibility  of  a 
weighted  combination  arises.  We  start  by  computing  what  has  been  found  to  be  a 
suitable  empirical  criterion,  given  by 

W  =  xQ)  +  1 . 7/(2  -  e/tt)  ,  (15) 

in  which,  since  x^^  <0  ,  the  two  terms  are  of  opposite  sign.  If  W  >  0  ,  we 
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use  Xqj  without  adulteration,  but  if  W  <  0  we  use  a  weighted  combination  of 
Xqj  with  Xq2  i  3°  that  the  starter  is  potentially  of  the  form 


II 

rO 

O 

*01 

1 

if 

W  i  0  j 

:03 

*01 

+  w(xQ2  -  X^) 

if 

W  <  0  ) 

-ally) 

the 

weight  w  is 

(-W)1/16 

,  the  16th  root  being 

(16) 


computed  as  /7YT  ;  the  transition  from  the  pure  x^^  to  the  weighted  is 

clearly  continuous. 


The  second  limitation  of  the  simple  starter,  (12),  was  associated  with 
(true)  values  of  x  close  to  -1,  the  potential  starter,  Xq^  >  giving  a  value 
much  too  close  to  -1,  This  flaw  has  been  dealt  with  by  applying  an  empirical 
factor  A  to  xQ3  ,  so  that  the  final  starter  (for  x  <  0  )  is  given  by 


x0  *  Ax03 


(17) 


The  formula  used  for  A  is 

A  -  1  +  ctx03(1  +  xQ1)  -  c2x03*/^1  +  W  •  (,8) 

with  values  of  and  equal  'empirically)  to  0.5  and  0.03  respectively! 

to  minimize  rounding  error,  1  +  in  (18)  is  computed  as  4/(4  +  T  -  T  )  , 

in  conformity  with  (12). 

A  point  concerning  the  patching  of  (12)  in  the  vicinity  of  (x,q)  3  (0,-1) 
is  worth  discussing.  We  sec  from  (15)  that  W  is  zero  at  this  point,  so  the 
patching  associated  with  (16)  can  have  little  effect  in  its  vicinity!  moreover, 
the  effect  of  the  patching  associated  with  (17)  is  very  slight  for  x  «  0  ,  so 
(12)  is  essentially  unpatched.  But  the  Xq^/Xq-j  patching  was  only  introduced 
because  xQ1  has  the  'wrong'  left-hand  derivative  at  (0,-1),  so  there  is  the 
appearance  of  a  contradiction  here.  The  paradox  is  resolved  if  we  bear  two 
points  in  mind.  First,  the  basing  of  (12)  on  the  value  of  Tq  means  that  there 
can  never  be  a  problem  when  x  »  0  .  Secondly,  so  long  as  q  is  not  exactly 
-1,  the  derivative  actually  has  the  'right'  value  (-4)  at  (0 , q) »  no  matter  how 
close  0  is  to  2m  .  Thus,  for  a  value  of  9  such  as  1,99999m,  the  unpatched 
starter  fails  in  the  vicinity  of  an  x-value  around  -0,05;  much  closer  to  zero, 
however,  all  is  well  again.  In  this  context,  the  inflexion-point  curve  of 
Lancaster  and  Blanchard  (Fig  4  of  Ref  10)  is  relevant,  though  their  associated 
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statement,  about  the  need  to  abandon  the  Newton-Raphson  process  in  the  region  in 
question,  is  unjustified  (as  already  noted  in  section  4), 

5.3  Multirevolution  starters 

For  m  >  0  ,  we  again  make  use  of  the  value  of  (corresponding  to 

x  ■  0  ,  for  the  given  q  ) ,  but  this  is  not  now  the  primary  source  of  information 

in  devising  starting  formulae.  The  primary  source  is  T.,  ,  the  minimum  value 

M 

of  T  for  the  given  q  ,  together  with  the  associated  value  of  x  ,  x^  .  As 

seen  from  Fig  2  (and  from  consideration  of  the  behaviour  of  T'  and  T"  ),  T„ 

and  arc  uniquely  defined;  further,  x^  ■  0  when  q  ■  1  ,  and  otherwise 

x^  >  0  .  Once  TM  is  known,  we  immediately  know  whether  the  Lambert  problem 

will  have  two  solutions  (T  >  T.,)  .  one  (T  =  T..)  ,  or  none  at  all  (T  <  Tu)  ,  but 

MM  M 

the  evaluation  of  (and  x^  )  is  a  non-trivial  matter,  itself  requiring  an 

iterative  process  -  a  solution  is  required  to  the  equation  T1  =  0  ,  For  this 
subsidiary  problem,  the  same  iteration  process  has  been  adopted  as  for  the 
Lambert  problem  (with  the  resulting  need  for  T'"  ,  as  remarked  in  section  4), 
so  there  are  two  topics  to  be  covered  here:  first,  the  starter  for  this  subsidiary 
problem;  and  secondly,  the  formulae  for  the  two  Lambert-problem  starters  that  will 
be  needed  as  soon  as  values  for  TM  and  xM  are  available,  assuming  that 
T  <  T  . 

In  regard  to  the  starter  for  x^  ,  we  need  the  value  of  ,  which  can  be 
recovered  from  q  ,  since,  from  (3)  and  (8), 

J0r  =  arg(2q  ,  1  -  q")  ,  (19) 

where  arg(x,y)  is  tan  * (y/x)  computed  unambiguously  and  implemented  in 
Fortran  by  the  ATAN2  function.  It  can  be  shown  that  for  -  it  (ie  for  q  “  0) 
a  good  approximation  to  x^  is  given  by  4/)3u(2m  +  1)}  ,  which  we  denote  by 
.  From  this  it  has  been  found,  empirically,  that  for  <  it  a  good 
starting  formula  is  (omitting  the  formal  zero-suffix  from  x^  ) 

XM  3  xM,1r(0r/’,)i  ’  C20) 

whilst  for  0  >  ir  it  is  (symmetrically) 

XM  “  XM,J2  • 


(21) 
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Having  determined  x^  and  the  corresponding  ,  we  now  suppose  that 

<  T  ,  so  chat  the  Lambert  problem  has  a  pair  o£  solutions  with  a  starting 

formula  required  for  both.  For  one  solution  we  have  x..  <  x  ,  so  x  is  certainly 

M 

positive,  and  this  leads  to  the  simpler  of  the  two  starters,  since  Tq  is  not 

used.  In  principle  we  get  x^  (the  required  starter)  by  inverting  the  particular 

bilinear  approximation  to  T  ,  as  a  function  of  (x  -  x^)2,  that  (i)  gives 

T  a  when  x  s  x^  ,  (ii)  makes  T  tend  to  infinity  as  x  tends  to  unity, 

and  (iii)  has  the  correct  second  derivative  (T")  at  x  *>  x„  .  This  bilinear  is 

M  M 

T  =  ’I’m  +  iT^(x  -  xM)2/!'  -  (x  -  xm)2/(I  -  xm)2|  ,  (22) 

ond  its  inversion  leads  to  the  starter 

*0  "  XM  +  [(T  ”  V^M  +  (T  ‘  rM)/(1  "  XM)ZJ]  ‘  (23) 

We  should  really  not  write  xQ  in  equation  (23),  as  it  requires  a  patch  that  is 
based  on  the  one  defined  by  (IV)  and  (18)  but  is  a  little  more  complicated;  the 
details  of  this  patch  are  omitted,  however,  but  Lhey  are  available  from  the 
listing  of  Appendix  D.  (A  new  constant,  c^  (  ,  is  involved,  as  well  as  the 

constants,  Cj  and  c2  ,  from  (18);  also,  m  and  0  are  arguments  of  the 

patch. ) 

For  the  starter  for  the  other  solution  we  have  x  <  x,,  ,  so  we  use  the  sign 

M 

of  T  -  to  distinguish  between  two  possibilities  (ignoring  the  third  possi¬ 
bility,  T  =  Tq  ,  which  is  as  trivial  as  the  one-solution  case,  T  “  ) .  If 

T  >  Tq  ,  then  x  <  0  and  we  proceed  exactly  as  when  m  »  0  and  x  <  0  ;  thus 
we  just  use  Tq  ,  not  requiring  T^  ,  and  we  patch  the  elementary  formula  for 
xQ  in  two  different  ways,  as  before  (with  the  second  involving  a  new  constant, 

^  ,  as  well  as  and  from  (18)).  Finally,  if  T  <  Tq  ,  we  have 

0  <  x  <  x..  ,  and  now  we  use  both  T  and  T„  .  We  base  x„  on  the  particular 
M  o  M  2. 

bilinear  approximation  to  T  ,  as  a  function  of  (x  -  x^)  '  ,  that  (i)  passes 

through  the  points  with  (x,  T)  equal  to  (x^  ,  T^)  and  (0  ,  T  ) ,  and  (ii)  has 

the  cor; act  second  derivative  (T")  at  the  former  point.  (We  could  match  T'  , 

M  o 

instead  of  ,  but  the  resulting  formula  would  be  as  complicated  as  equation 

(25)  following.)  This  bilinear  approximation  is 


*T"(x 

M 


V 


I  + 


(*  -  vw<t0  -  y 


(24) 
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and  its  inversion  leads  to  the  starter 


(  t-’h 

iTM  “  (T  ’  TM){lTM/To 

-  v  “ 

(25) 


equation  (25)  is  a  complicated  formula,  but  no  patching  of  this  starter  has  been 
found  necessary. 

6  COMPLETION  OF  CONVERGENCE 

29-3  j 

It  was  decided  (again  following  the  work  on  the  two  Kepler  equations) 

to  aim  at  a  fixed  number  of  iterations  for  x  ,  rather  than  employ  a  convergence 
test.  In  program  testing,  carried  out  on  a  PRIME  computer  providing  14  decimal 
digits  in  double  precision,  it  was  found,  for  m  *  0  ,  that  three  iterations 
always  sufficed  for  the  determination  of  x  to  13  digits,  suggesting  that  (for 
levels  of  precision  up  to  this)  no  truncation  error  remained;  it  was  legitimate, 
therefore,  to  fix  the  number  of  iterations  at  three.  The  foregoing  is  an  over¬ 
simplified  statement  of  the  accuracy  achieved,  however,  and  an  amplified  version 
now  follows. 

It  was  arranged,  in  testing,  that  the  relative  error  in  x  (true  value 
assumed  known  -  see  section  8),  after  a  pre-set  number  of  iterations,  should  be 
computed,  as  well  as  the  relative  error  (residual)  in  the  value  of  T  finally 
computed  by  TI.AMB  (as  compared  with  the  given  T  ) ;  the  smaller  of  these  two 
relative  errors,  e  say,  was  registered  for  each  test  case,  following  the 

rationale  given  in  section  1  of  Ref  29,  and  after  three  iterations  it  was  found 

-13 

that  c  never  exceeded  10  .  (This  rationale,  to  be  invoked  again  in  section  8, 

is  based  on  the  proposition  that  the  accuracy  in  a  numerical  solution  of  the 
general  equation  f(x)  =»  Y  should  always  be  assessed  in  terms  of  the  numerically 
smaller  of  the  relative  error  in  x  and  the  relative  residual  in  Y  ;  assessment 
in  terms  of  just  the  former  could  amount  to  a  demand  for  the  impossible.)  For 
completeness,  it  is  worth  remarking  that  when  the  process  was  reduced  to  two 
iterations  the  maximum  value  of  e  was  found  to  be  about  2  x  10  this  being 
the  relative  error  in  x  that  arises  when  0  w  1 . 7 1  it  and  x  «  -0.9944  ,  For 

-3 

just  a  single  iteration,  on  the  other  hand,  the  maximum  value  of  e  is  4.3  *  10  , 
occurring  for  values  of  0  approaching  2x  and  for  x«  2.28.  Finally,  the 
maximum  value  of  e  prior  to  any  iteration,  ie  due  to  the  starter  itself,  is 
about  0.5,  being  associated  with  the  100%  over-estimate  by  xq  ,  when  0  is 
small,  noted  in  section  4.  Consideration  of  the  convergence  when  m  >  0  is  held 
over  to  section  8. 
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For  ths  convergence  of  x^  in  the  'subsidiary  problem'  (as  defined  in 

section  5.3),  it  was  soon  clear  that  a  fixed  number  of  iterations  would  not  be 

right.  There  was  something  of  a  dilemma  here,  in  regard  to  the  philosophy  to 

adopt,  since,  on  the  one  hand,  great  accuracy  in  should  not  really  be 

required,  as  x^  is  only  a  step  to  the  starter  for  the  main  problem;  in  critical 

cases,  on  the  other  hand,  an  incorrect  could  lead  to  the  erroneous  con- 

elusion  that  the  particular  Lambert  problem  possesses  no  solution.  It  was 

eventually  decided  to  op  .  e  on  the  (conservative)  basis  that  iteration  would 

end  as  soon  as  the  value  of  x^  changed  by  less  than  throe  parts  in  10  during 

the  current  iteration.  With  the  cubic  convergence  of  HaLley's  method,  this  meant, 

in  principle,  that,  in  a  hypothetical  further  iteration,  x  would  not  change  by 

20  ” 

more  than  three  parts  in  10  -  this  follows  from  r.he  discussion  of  convergence, 

in  the  context  of  Kepler's  equation,  in  Ref  31,  l‘hi3  expectation  was  confirmed 
when  it  was  found  that  there  was  usually  no  change  in  the  value  of  T  itself  if 
the  'further  iteration'  was  actually  performed.  Tests  for  values  of  0  covering 
the  full  range  of  q  were  carried  out;  the  maximum  number  of  iterations  needed 
to  satisfy  the  x^  criterion  was  found  to  be  nine,  but  the  number  was  only  three 
in  the  vast  majority  of  cases,  This  was  true  even  for  largo  values  of  0  ,  the 
convergence  being  essentially  dependent  on  q  rather  than  m  .  To  provide  a 
guaranteed  exit,  however,  it  was  arranged  that  if  the  criterion  was  not  satisfied 
within  twelve  iterations,  then  the  process  would  be  abandoned  and  a  flag  set; 
reference  to  this  has  already  been  made  (section  4)  . 

It  is  worth  noting  how  convergence  for  x^  would  be  affected  by  substitut¬ 
ing  the  Newton-Raphson  process  for  the  Halley  process.  For  most  of  the  q-range, 
the  effect  would  be  to  add  only  a  single  Iteration,  taking  the  total  from  three 
to  four.  For  values  of  q  approaching  1,  however,  ie  for  correspondingly  small 
values  of  ,  there  is  a  steady  rise  in  the  number  of  iterations  required  by 
the  Newton-Raphson  process  -  the  maximum  number  experienced  in  the  PRIME  testing 
was  16,  corresponding  to  nine  for  the  Halley  process.  To  avoid  this  behaviour, 
it  would  be  necessary  to  improve  the  starting  formula,  (20), 

7  COMPUTATION  OF  VELOCITY 

We  come  at  last  to  the  solution  of  an  individual  Lambert  problem.  It  is 
assumed  that  we  have  obtained  a  value  of  x  (and  also,  where  appropriate,  x+  ) 
for  the  Lambert-equivalent  class  of  problems,  and  it  remains  to  compute,  for 
both  of  the  points  P.  and  P,  ,  the  velocity  components  (VD  and  V„)  correspond- 
ing  to  the  specific  quantities  (r^,  r^,  9  and  At)  of  the  original  problem.  This 
computation  is  performed  by  the  overall  I.ambert-solving  subroutine,  VLAMB,  which 
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calls  XLAMB  (which  itself  calls  TLAMB) .  Vl.AMB  is  listed  in  Appendix  Ej  it  will 

be  seen  that  its  arguments  are  precisely  the  quantities  summarized  (for  the  full 

Lambert  solution  procedure)  in  section  1.  (For  'minimalist'  derivation  of  the 

formulae  for  VD  and  V_  ,  see  Ref  28.) 

K  i 

A  formula  for  V_  ,  was  given  by  Lancaster  and  Blanchard'^,  This  can 
K|  * 

be  improved  upon,  however,  to  make  it  more  accurate  and  efficient.  The  resulting 
formula  is 

Vr  !  =  v|(qz  -  x)  -  o(qz  +  x)|/r1  ,  (26) 

where  y  -  /(ps/2)  and  p  =  (rt  -  r^l/c  .  The  formula  for  Vr  ^  >  similarly,  is 

Vr  2  *  "  v|(qz  -  x)  +  p (qz  +  x)}/r2  •  (27) 

28 

Sarneeki's  interpretations  of  x  and  z  are  immediately  apparent,  on  putting 

0  =  1. 

We  must  have  | o I  <  1  ,  of  course,  except  that  p  is  indeterminate  when 

c  *  3  (P1  and  coincide),  in  which  case  VLAMB  arbitrarily  sets  it  to  zero. 

But  if  c  a  0  ,  q  *  ±1  ,  so  that  z  =  |x|  by  (7),  and  then  for  x  )  0  there  are 

two  possibilities,  depending  on  whether  q  and  x  are  of  the  same  or  opposite 

sign.  Only  when  the  signs  are  the  same  do  V  ,  and  Vu  ,  themselves  become 

indeterminate  (see  also  the  explanatory  material  in  sections  1,  2.1  and  3),  and 

then  the  arbitrary  value  of  c  yields  conventional  values  for  V„  .  and  V  „  . 

R)  I  R j  z 

There  is  no  indeterminacy  when  q  and  x  have  opposite  signs,  on  the  other  hand, 
since  VD  .  and  VB  ,  (also  V„  .  and  V_  B  )  are  then  independent  of  p  j  we 

K#  I  K|  ^  I  f  1  1 1  ^ 

now  have  the  rectilinear  orbit3  illustrated  by  (the  limiting  form  of)  Fig  4a  (on 
the  right)  and  Fig  4e  (on  the  left). 

A  formula  for  V„  was  not  given  directly  in  Ref  10,  but  only  via  the 
*■  »  * 

orbital  elements  a  and  e  ,  with  a  serious  threat  to  accuracy  in  awkward  cases. 
The  direct  formula  is 

VT  (  “  ya(z  +  qx)/r(  ,  (28) 

where  in  principle  a  is  defined  as  (1  -  p2)  ^  i  to  minimize  rounding  error, 
however,  in  practice  we  compute 
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setting  a  a  i  (compatible  with  p  )  when  c  “  0  .  Similarly, 

Vj  j  "  yo(z  +  qx)/r2  .  (30) 

Naturally,  V^,  ^  and  2  a-e  indeterminate  in  the  same  circumstances  as 
j  and  VR  2  ,  with  compatible  resolution  of  this  indeterminacy. 

for  use  in  equations  (26)  to  (28)  and  (30),  the  values  of  qz  -  x  ,  qz  +  x 
and  z  +  qx  are  given,  to  optimum  accuracy,  as  special-case  output  from  the 
subroutine  TLAMB ,  in  a  post-XLAMfl  direct  call  by  Vl.AMB.  The  special  case  is 
signalled  by  the  setting  of  (input)  N  to  -I,  which  causes  the  quantities  in 
question  to  be  computed  in  place  of  T'  ,  T"  and  T"'  (see  Appendices  li  and  C)  . 

8  TESTING  RATIONALE  AND  RESULTS 

The  subroutines  TLAMB,  XLAMB  and  VLAMB  were  tested  for  a  wide  range  of 
data,  starting  with  TLAMB  which  is  a  self-contained  procedure  implementing  a 
direct  algorithm  (for  computing  T  from  x  ).  Problems  of  accuracy  with  this 
subroutine  seemed  most  likely  to  arise  with  input  for  which  q  was  close  to  ±1, 
or  else  x  was  close  to  ±1  or  zero,  so  testing  was  particularly  thorough  for 
data  in  these  categories.  Further,  there  are  (for  m  »  0)  two  transitions  (from 
ellipse  through  'series'  to  hyperbola),  one  each  side  of  x  =  1  (see  Appendix  B) ; 
the  regions  in  the  vicinity  of  these  transitions,  defined  (in  TLA11B)  in  practice 
by  x  =  A). 6  and  x  =  A 74  ,  were  tested  the  most  carefully  of  all,  with  con¬ 
sistency  carefully  monitored.  The  tests  were  entirely  satisfactory,  and  Fig  2 
was  based  on  the  output  -  its  consistency  with  the  Figure  given  by  Lancaster 
at  al  was  an  additional  confirmation  that  all  was  in  order. 

The  testing  of  XLAMB  followed  a  natural  procedure  for  validating  the 
iterative  solution  of  an  equation.  Thus,  the  parameters  of  this  subroutine  are 
effectively  0  and  T  ,  but  instead  of  ranging  over  T  the  test  data  actually 
ranged  over  x  .  Then  each  'true'  x  (with  a  given  0  )  was  the  source  of  a 
nominal  test  value  of  T  (via  TLAMB),  after  which  the  testing  of  XLAMB  could 
proceed,  with  the  true  x  used  (at  the  end  of  the  test)  merely  as  an  accuracy 
evaluator.  Only  one  solution  (in  two-solution  cases)  could  be  relevant  to  this 
evaluation;  further,  the  specific  one-solution  cases  and  the  no-solution  cases 
then  required  separate  testing.  The  general  testing  was  very  thorough,  with 
values  of  0^  taken  very  close  (and  even  equal)  to  0  and  2ir  ,  and  (as  already 
remarked  in  section  6)  it  has  indicated  that,  for  m  “  0  ,  13-digit  accuracy  is 
always  achieved  within  three  iterations.  (The  testing  extended  to  a  value  of 
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10  for  x  ;  this  confirmed  that  full  accuracy  was  maintained  as  the  attractive 
force  was  effectively  reduced  to  zero.) 

A  particular  case  was  examined  in  detail,  namely,  that  for  0  ~  1.9936tt 

and  x  »:  0.159  ,  chosen  because  this  implied  an  input  of  q  ■  0.99  and  T  ■  5.0, 

.  23 

for  which  the  paper  by  Battin  and  Vaughan  indicated  the  worst  convergence 

behaviour  of  the  cases  the  authors  considered.  In  fact  Ref  23  registers  14 

iterations  as  necessary  in  the  modified  version  of  Gauss's  method  of  alternating 

substitution  (see  Appendix  A),  just  to  get  a  solution  to  8-digit  accuracy. 

24 

(Boltz  claims  to  achieve  faster  convergence  Chan  Battin  and  Vaughan,  but  be 

does  no  better  in  the  worst  cases  compared.)  Ref  23  suggests  that  methods  of  the 

Newton-Raphson  type  fail  in  these  circumstances,  but  this  is  far  from  so  (in 

spite  of  the  remark  of  Lancaster  and  Blanchard'^  that  originated  this  suggestion, 

already  commented  upon  here  in  sections  4  and  5.2).  Thus  the  values  of  e 

(from  section  6)  after  zero,  one  and  two  iterations  are,  respectively,  0.043, 

-4  - 1 2 

1.3  x  10  and  5.2  *  10  ;  after  three  iterations,  c  is  too  small  to  detect, 

-34 

but  it  would  evidently  be  around  4  *  10  in  the  absence  of  rounding  error. 

The  figures  quoted  are  for  the  Halley  iteration  process  incorporated  in  XLAMB, 

but  (as  it  happens)  the  results  are  no  worse  if  the  Halley  process  is  replaced 

hy  the  Newton-Raphson  process.  However,  this  conspicuously  good  behaviour  is 

largely  due  to  the  bilinear  starter,  so  it  was  worth  seeing  what  would  happen  if, 

for  compatibility  with  Ref  23,  we  took  xn  from  its  value  for  the  circular 

U  7-1 

orbit  through  and  P,,  .  The  formula  for  this  is  -  q(1  +  q")  , 

giving  about  0.704  for  the  case  considered.  Then  the  values  of  e  ,  with  the 

“6 

built-in  Halley  process,  increased  to  0.36,  0.012  and  3,4  x  10  ,  for  iterations 

up  to  the  second,  the  value  after  three  iterations  still  being  '  submerged  in 
rounding  noise  '.  When  the  Halley  process  was  replaced  by  the  Newton-Raphson, 

however,  it  was  another  story:  values  of  c  ,  taken  now  as  far  as  the  third 

-4 

iteration,  were  0.36,  0.23,  0.034  and  6,2  x  10  .  For  this  example,  the  con¬ 
clusion  is  clear:  only  if  the  starting  formula  and  the  iteration  process  are 
both  degraded  will  convergence  deteriorate  seriously:  even  then  our  normal 
criterion  would  be  met  after  at  most  another  two  iterations. 

-13 

For  m  >  0  ,  the  value  of  c  is  no  longer  less  than  10  in  all  cases, 

-13 

though  for  m  =  1  the  maximum  value  obtained  was  still  only  1.1  x  10  .  For 

m  >  1  ,  e(max)  grows  steadily  -  the  initial  growth  rate  of  its  common 

logarithm,  with  respect  to  m  ,  is  about  0.3,  but  the  rate  falls  off  gradually, 
c(max)  being  about  3  *  12  ^  for  »  =  5  ,  8  *  10  11  for  ni  =  10  ,  6  x  10  ®  for 
m  =  30  and  3  x  io  5  for  m  =  100  .  (Up  to  about  this  point,  full  accuracy  would 


26 


be  restored  by  going  to  a  fourth  iteration,  since  this  reduces  e(max)  to  about 

1,7  x  10  when  m  =  100  ;  it  rises  to  about  2.4  »  10  ^  when  m  =  200  ,  however,) 

The  growth  of  e(max)  with  ra  is  entirely  due  to  a  defect  in  the  patching,  for 

x  <  0  ,  such  that  the  condition  for  invocation  of  the  patch  involving  2  is 

sometimes  not  satisfied  when  one  would  like  it  to  be.  As  the  patching  is 

empirical,  it  could  in  principle  be  improved  without  much  difficulty,  but  the 

need  to  do  this  is  so  obviously  academic  that  the  matter  has  not  been  pursued. 

It  is  worth  remarking,  however,  that  for  large  values  of  m  ,  the  accuracy  of 

XLAMB  is  surprisingly  sensitive  to  the  particular  values  assigned  to  c^  and 

c.  small  changes  can  easily  lead  to  a  worsening  in  accuracy  by  several 

^  »  * 

orders  of  magnitude.  In  view  of  the  empirical  nature  of  the  patching,  it  is 
perhaps  to  be  wondered  at  that  the  multirevolution  starters  have  performed  as 
well  as  they  have. 

The  testing  of  VLAMB  was  based  on  the  idea  that,  for  each  Lambert  solution, 
the  output  V  ,  and  V  .  of  the  subroutine  can  be  combined  with  r  (together 
with  a  polar  reference  angle,  taken  as  zero)  to  constitute  the  four  necessary 

components  of  data,  associated  with  the  point  P  ,  for  input  to  the  subroutine 

27  ' 

PV2ELS.  This  subroutine  then  generates  the  corresponding  set  of  four  universal 

two-dimensional  orbital  elements,  one  of  which  is  t  ,  the  time  from  pericentre 

(a  conventional  paint  if  the  orbit  is  circular).  When  t  is  updated  to  t  +  At  , 

with  the  other  three  elements  unchanged,  a  nominal  position  (plus  velocity)  for 

the  point  can  be  obtained  from  ELS2PV,  the  subroutine  that  is  inverse  to 

PV2ELS.  This  position  is  specified  as  a  radius  vector  and  a  polar  angle,  so  the 

performance  of  VLAMB,  and  hence  the  success  of  the  overall  solution  procedure, 

can  be  assessed  by  a  direct  comparison  of  these  quantities  with  the  input  (test) 

values  r2  and  6  j  we  denote  the  (absolute  values  of)  the  differences,  in  r2 

and  0  ,  by  6r  and  60  .  Ideally,  looking  at  relative  errors,  we  should  like 

6r/r„  and  69/0  ,  which  we  can  denote  by  6  and  6.  ,  to  be  no  greater  than, 
l  _ ,  *3  r  u 

say,  5  *  10  (on  the  PRIME  computer). 

For  various  reasons,  this  'ideal'  requirement  (in  VLAMB  testing)  is  too 
stringent.  First,  if  the  velocity  is  (in  relative  terms)  very  great  at  , 

then  a  Large  magnify cation  of  the  relative  rounding  error  may  occur  that  is 
completely  unavoidable.  For  6^  ,  this  magnification  will  be  allowed  for,  in 
principle,  if  we  do  not  automatically  divide  6t  by  ,  but  instead  use  rR  , 
where 

rR  -  max(r2  ,  |VR  2|  At)  .  (31) 
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(We  define  r  similarly,  and  apply  the  factor  r^/r^  to  56  ■)  This  does  not 
fully  compensate  for  the  velocity  effect,  however,  since  At  may  have  a  large 
derivative  with  respect  to  the  value  of  x  determined  by  XLAMBi  this  secondary 
effect  is  allowed  for  if  we  replace  At  in  (31)  by  max(At  ,  |x  d(At)/dx]  ), 

(The  rationale  here  is  the  same  as  the  one  involved  in  the  definition  of  e  in 
section  6;  the  value  of  the  derivative  comes  at  once  from  T'  ,  using  the 
relation  between  At  and  T  in  equation  (1),) 

With  the  foregoing  allowance  for  the  velocity  at  P2  ,  the  modified  value 
of  6  satisfied  the  'ideal  requirement'  for  all  tests  with  r^  ,  but  a 

further  complication  with  the  test  criteria  was  observed  in  tests  with  r,  >  . 

(We  have  hitherto  thought  of  as  being  at  least  as  great  as  r?  ,  because  of 

the  association  with  the  definition  of  x  ,  but  this  restriction  did  not  apply 
in  the  testing  of  VLAMB.)  The  difficulty  may  be  understood  in  relation  to  situ¬ 
ations  with  r2  ">:>  ri  >  au<^  a  value  of  x  corresponding  to  a  semi-major  axis 

with  a  >>  ,  since  then  the  determination  of  the  element  a(»p/a)  by  PV2ELS 

is  inevitably  inaccurate.  In  principle,  it  would  be  better  to  replace  it  by  the 
value  given  (from  x  and  s  )  by  equation  (1) ,  but  as  the  object  of  the  exercise 
was  to  test  the  new  subroutine  (VLAMB)  with  the  existing  subroutines  PV2ELS  and 
ELS2PV,  this  has  not  been  done.  Instead,  the  question  "what  empirical  further 
relaxation  in  the  definition  of  5^  (and  similarly  60  )  would  cause  the  diffi¬ 
culty  to  disappear?"  was  addressed.  It  was  found  that  an  additional  division, 
in  computing  the  relative  errors,  by  the  ratio  ^/r^  '  or  rather  by 
max(1,  r2/r|)  t0  cover  the  other  case  (r(  >,  r?)  as  well,  would  for  the  most 

part  resolve  the  problem;  complete  resolution  became  possible  on  replacing  the 

1  3 

empirical  ratio  xj^X\  *'y  the  even  more  empirical  (r^/r^  . 

With  5  and  5  adjusted  on  the  rationale  of  the  last  two  paragraphs,  it 

r  0 

was  found  that  5  always  satisfied  the  'ideal  requirement1,  but  that 
could  fail  to  do  so  when  0  was  less  than  about  1  radian.  Under  these  circum¬ 
stances,  replacement  of  the  relative  error  (50)  by  the  absolute  error  (60) 
permitted  the  requirement  to  be  satisfied,  but  this  is  hardly  surprising,  and, 
in  v iew  of  the  desirability  of  maintaining  relative  accuracy  for  Lambert  problems 
with  small  values  of  3  and  At  ,  it  was  important  to  know  whether  the  loss  of 
accuracy  arose  ir.  the  Lambert  procedure  itself  or  in  the  test  procedure  via 
PV2ELS  and  ELS2PV.  Tests  for  very  small  values  of  both  the  input  parameters 
have  indicated  that  the  errors  are  entirely  in  the  test  procedure.  The  expla¬ 
nation  is  that  PV2ELS  and  ELS2PV  were  not  designed  for  extreme  accuracy  in  moving 
between  close  points  on  a  given  orbital  path;  in  particular,  this  follows  from 
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the  subroutines'  built-in  reference  to  a  purieentre,  which  is  likely  to  be  a 
point  that  is  remote  from  the  given  points  that  are  neighbours,  (it  would  bo 
possible  to  formulate  the  subroutines  in  a  different  way,  but  this  would  involve 
considerable  complication  in  the  solution  of  the  two  forms  of  Kepler's  equation; 
thus,  for  an  elliptic  orbit  it  would  be  necessary  to  solve  an  equation  for  the 

eccentric-anomaly  difference  E.  -  E,  ,  rather  than  soLving  the  classical,  equation 

1 13  i 

for  E^  .  Lancaster  has  written  on  the  Lambert  problem  for  short  arcs,  whilst 
Battin'2  has  considi red  the  effective  inversion  of  Lambert's  boundary-value 
problem  to  provide  a  solution  for  Kepler's  initial-value  problem!) 

With  the  original  definitions  of  6  and  6^  modified  as  indicated, 

r  0 

including  the  replacement  of  the  latter  by  60  when  0  is  less  than  I,  their 
.  -  1 3  . 

values  remained  below  5  x  10  in  all  the  tests  carried  out.  In  this  testing, 
the  same  range  of  values  for  0  was  used  as  in  testing  VLAMB,  with  the  same 
range  of  input-x  used  as  source  for  the  nominal  values  of  At  ,  Separate  values 
of  r(  and  r 2  were  now  provided,  though  only  their  ratio  was  actually  signifi¬ 
cant;  this  meant,  of  course,  that  the  input  0  was  no  longer  identical  with  0  . 

Various  values  for  r„/r.  were  selected,  covering  the  range  from  10  ^  to  106, 

V  .  '  - 1 3 

so  the  fact  that  (modified)  6^  and  6^  could  be  held  to  3  x  10  ,  over  such 

a  wide  range  of  x  ,  0  and  r(/r2  '  must  regarded  as  an  entirely  successful 

outcome  to  the  VLAMB  testing;  it  incidentally  increased  the  confidence  placed  in 

the  robustness  of  PV2ELS  and  ELS2PV,  except  in  regard  to  the  point  of  the 

preceding  paragraph. 

9  CONCLUSIONS 

A  study  of  the  literature  on  Lambert's  orbital  boundary-value  problem  shows 
that  the  crucial  contribution  to  the  subject  was  made  by  Lancaster,  Blanchard 

g 

and  Devaney  in  1966.  Some  of  the  recent  papers,  unfortunately,  do  not  refer  to 
Ref  8  at  all,  and  their  treatment  of  the  problem  is,  to  this  extent, 
retrogressive. 

8 

The  original  paper  of  Lancaster  ct  al  was  very  short  (less  than  ll  pages) 
and,  though  their  approach  was  amplified  in  a  later  paper'^,  there  has  been  a 
need  for  the  approach  to  be  extended  to  a  general  computing  procedure  for  solving 
the  Lambert  problem.  The  present  Report  fills  this  gap,  and  has  addressed,  in 
particular,  the  following  topics  not  covered  by  Lancaster  at  all  the  provision 
of  universal  starting  formulae  for  the  procedure's  iteration  process;  the  method 
of  iteration;  the  minimization  of  rounding  error  in  all  circumstances  (not  just 
around  the  parabolic-transition  region);  the  accurate  computation  of  velocity; 
and  the  test-validatior  of  the  procedure  (with  emphasis  on  extreme  cases). 
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The  overall  conclusion  is  that  the  procedure  provides  excellent  accuracy, 
in  an  efficient  manner,  for  any  problem  that  might  arise  in  practice,  Three 
iterations  of  the  iteration  process  suffice  for  the  effective  elimination  of 
rounding  error,  when  working  to  an  accuracy  of  not  more  than  )3  significant 
figures.  When  only  5-  or  6-figure  accuracy  is  necessary,  two  iterations  will 
be  adequate  -  in  most  cases  they  will  be  a  great  deal  more  than  merely  'adequate'. 

The  Report  has  attempted  to  facilitate  a  deeper  understanding  of  Lambert's 
problem  by  introducing  the  concepts  of  L-similarity  and  L-congruence  (for  the 
basic  triangle  involved  in  the  problem)  and  of  Lambert  invariance  (for  physical 
parameters) . 
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Appendix  A 

GAUSS'S  METHOD  AND  XHE  PAPER  BY  BATTIN  AND  VAUGHAN 

Ref  23  represents  the  culmination  of  a  long  study  of  Lambert's  problem  by 

Battin  and  various  co-workers.  Ref  2  (his  text-book)  and  Refs  12,  18,  19  and  23 

do  not  exhaust  the  list  of  Battin's  contributions  to  the  subject,  but  Ref  23 

represents  a  change  of  approach  in  which  the  authors  return  to  the  classical 

32  33 

method  of  Gauss.  The  essence  of  Gauss's  method,  of  which  Moulton  and  Plummer 
give  traditional  text-book  accounts,  lies  in  its  iterative  determination,  by 
'alternating  substitution'  (not  a  standard  term),  of  a  pair  of  Lambert-invariant 
parameters,  traditionally  denoted  by  x  and  y  .  For  an  elliptic  orbit,  Gauss's 
x  is  sin  ir  (E^  “  E^)  ,  where  Ej  and  E,  are  the  eccentric  anomalies  of  the 
points  Fj  and  P,,  j  x  also  identifies  with  4(1  —  g)  ,  where  g  is  defined 
by  equation  (B-9)  of  Appendix  B.  Again,  y  is  a  quantity  defined  by  the  Lambert 
triangle  CP(P2  an^  t*ie  re(lu ir ed  orbital  path,  being  the  ratio  of  the  area  of 
the  'curvilinear'  triangle  to  that  of  the  ordinary  'linear'  triangle,  where  the 
curvilinear  triangle  has,  for  its  'side'  ,  the  orbital  arc  rather  than  the 

chord  (so  y  is' infinite  when  0  =  ir  )  i  this  identifies  Gauss's  y  with  T/4qa 
in  the  notation  oi  Appendix  B,  where  a  is  given  by  equation  (B-4) .  Two  equations 
are  available  for  connecting  x  and  y  ,  via 

y3  -  y2  =  mX  (A- 1 ) 

and 

x  -  m/y2  -  l  .  (A-2) 

In  ( A— 1 )  and  (A-2),  l  is  equivalent  to  q  (in  the  notation  of  the  main  text) 
since 

£.  =  (1  -  q)2/4q  (A-3a) 

2  • 

=  sin  tO^/cos  i  ©r  .  (A-3b) 

Computation  of  £  from  0  ,  via  (A-3b),  is  more  accurate  than  computation  from 
q  ,  via  (A-3a) ,  when  0^  is  small.  For  computation  from  ,  rather  than  0  , 

Gauss  recommended  a  formula  that  may  be  derived  from  the  basic  relation 

/(r^)  cos  J  0  =  qs  =  4 (r  1  +  r,,)  cos  j  0  , 

I 

i 
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which  is  just  a  combination  of  equations  (6a)  and  (6b)  of  the  main  text.  (This, 
incidentally,  identifies  Q  with  2f  in  Ref  17.)  Gauss's  formula  is 

2  i  5 

A  ■  (sin  r6r  +  tan  2ie)/cos  5  6  (A-3c) 


where  u)  is  given  by 

l 

tan(w  +  £rr)  ■  (r2/r1)ir  , 

and  it  is  evident  that  (A-3c)  is  just  a  generalization  of  (A-3b) .  However,  we 

,  ,  2 
can  get  a  still  more  accurate  value  of  l  if  we  replace  tan  2w  in  (A-3c)  by 

2 

(/r^  -  /r^)  /ft/Cr^)  ,  and  we  do  best  of  all,  when  an  accurate  value  of 
| f ^  —  r ^ I  is  available,  if  the  replacement  is  by 

( r ,  ~  r2)2/{4/(rlr2)(/ri  +  /?2>2}  • 

To  return  to  the  quantities  introduced  in  (A-1)  and  (A-2),  m  (actually  written 
m2  by  Moulton,  following  Gauss's  Theoria  Motus34)  is  closely  related  to  T 
(present  main  text)  since 


m  »  u(At)2/(8s3q3)  ,  (A-4a) 

from  which,  using  equation  (2)  of  the  main  text,  we  get 

m  “  T2/64q3  .  (A-4b) 

Finally,  X  is  a  function  of  x  alone,  most  concisely  expressed,  in  terms  of 
the  usual  hypergeometric  function,  as  s’F (3 , 1 , 2 i  ;x)  ,  but  most  effectively  com¬ 
puted  using  continued  fractions.  The  determination  of  y  from  x  in  (A-1)  then 
involves  the  solution  of  a  cubic  equation. 

Gauss  based  the  iterative  determination  of  x  and  y  ,  to  satisfy  (A-1) 
and  (A-2)  simultaneously,  on  a  method  of  successive  substitution  (a  form  of 
relaxation)  that  can  more  descriptively  be  referred  to  as  'alternating  substi¬ 
tution'  when  only  two  equations  are  involved.  To  see  how  the  method  operates, 
consider  (A-1)  and  (A-2)  generalized  to  y  =  f(x)  and  x  =  g(y)  respectively, 
with  an  initial  estimate  available  for  x  s  then  we  compute  y^  «  f(xQ)  * 
x2  "  »  y3  ”  £(Xj)  et.' ?,  and  need  to  assess  how  the  alternating  sequences 

converge,  if  at  all.  So  let  £'  and  g'  denote  the  derivatives  of  f  and  g  , 
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and  suppose  Chat  £'  and  g'  ate  reasonably  constant  over  the  square  defined 
by  (Xq  ,  y^  and  the  true  solution  (x  ,  y) .  Then  if  | f 1  g *  |  <  1  ,  an  error,  c 
aay,  in  x^  leads  to  an  error  of  order  f'g’e  in  x^  .  On  this  assumption,  the 
process  is  linearly  convergent!  the  rate  of  convergence  is  rapid  if  |fg|  <<  1  , 
but  still  only  linear.  Thus,  though  the  process  is  likely  to  be  very  robust,  it 
is  inherently  inferior  to  the  standard  iteration  processes  for  solving  a  single 
equation  in  a  9ingle  variable  -  in  particular  the  Newton-Raphson  process  which 
gives  quadratic  convergence  and  the  Halley  process  which  gives  cubic  convergence. 
The  main  weakness  in  Gauss's  method  comes  from  the  infinities  in  £  and  m  when 
q  ■  0  (which  should  bo  the  simplest  of  all  cases),  and  a  subsidiary  defect  lies 
with  the  slow  convergence  in  other  cases.  Battin  and  Vaughan  have  improved  the 
method  significantly  in  both  respects,  dealing  with  the  infinities  by  a  redefi¬ 
nition  of  all  the  quantities  £  ,  m  ,  x  ,  y  and  X  ,  in  such  a  way  that  (A-1) 
and  (A-2)  are  (formally)  still  the  equations  to  be  solved.  In  particular,  the 
redefinition  of  m  is  as 

m  **  T2/ ( 1  +  q)6  ( 

this  still  leaves  an  infinity,  but  now  it  is  for  q  ■*  -1  (0  a  multiple  of  2r  , 
in  the  notation  of  the  main  text),  which  is  a  lesser  fault  but  a  fault  nonetheless, 
Battin  and  Vaughan  deal  with  the  other  weakness  by  an  adjustment  of  the  equations 
(A-1)  and  (A-2)  themselves,  in  a  way  that  leads  to  a  dramatic  speeding  up  of  the 
convergence  (via  a  reduction  in  the  value  of  |f'g'|  ,  though  they  do  not  express 
the  alternating  substitution  process  in  this  way). 

It  may  well  be  supposed,  however,  that  a  much  more  drastic  modification  of 
Gauss's  method  could  be  devised.  If  we  remove  the  denominator  in  (A-5) ,  and  then 
formulate  (A-1)  and  (A-2)  as  an  equation  in  a  single  unknown  (as  is  hinted  at  in 
Ref  14),  we  remove  both  the  deficiencies  completely,  In  doing  this  in  the  most 
efficient  way,  however,  we  effectively  just  recover  the  method  of  Lancaster  et  at, 
and  it  is  ironic  that  the  only  reference  to  this  method  made  by  Battin  and  Vaughan 
is  a  remark  to  the  effect  that  "the  difficulties  arising  when  6  is  a  multiple  of 
2u  were  recognized  by  Lancaster  and  Blanchard1^  and  necessitate  the  abandonment 
of  the  Newton-Raphson  process".  It  is  perfectly  true  that  the  possibility  that 
this  process  might  have  to  be  abandoned  was  suggested  in  Ref  10.  In  reality, 
however,  as  should  be  clear  from  section  S  of  the  present  Report,  this  is  the 
last  circumstance  in  which  it  is  actually  profitable  to  change  the  iteration 
process;  the  Gauss  process  is  inevitably  still  at  its  worst,  in  spite  of  the 
improvements  of  Battin  and  Vaughan,  whereas  the  Newton-Raphson  process  (or  better 


027 


Appendix  A 


33 


the  Halley  process)  continues  to  do  extremely  well  in  three  iterations!  this 
success  must  be  attributed  to  the  use  of  a  good  starting  formula  for  the  x 
Lancaster  ct  at. 


of 
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Appendix  B 

DETAILS  _0F_ THE  TI.A1-IB  ALGORITHM 

The  normal  function  of  the  subroutine  TLAMB,  which  is  Listed  in  Appendix  C, 

is  to  compute  the  value  of  the  Lambert-invariant  parameter  T  ,  defined  in 

section  3  of  the  main  text,  together  with  as  many  of  its  first  three  derivatives 

(with  respect  to  the  Lambert-invariant  parameter  x  )  as  are  required.  The  input 

2 

arguments  of  the  subroutine  are  m  ,  q  ,  1  -  q  ,  x  and  n  ,  whore  n  ,  when 
0  <  n  <  3  (normal  operation),  specifies  the  number  of  derivatives  to  be  output, 
the  output  arguments  being  T  ,  T'  ,  T"  and  T"'  .  When  n  =  -1  ,  the  sub¬ 
routine  has  a  special  function,  nameLy,  to  compute  the  quantities  S  ,  B  and  A 
that  will  be  defined  in  due  course;  these  quantities  then  replace  the  arguments 
T'  ,  T"  and  T"1  ,  there  being  no  output  T  . 

The  basic  algorithm  for  T  is  unchanged  from  the  definitive  work  of 
8  10 

Lancaster  at  al  '  ,  and  we  use  the  same  notation,  except  for  the  introduction 

of  a  ,  A  ,  6  and  B  ,  and  for  the  use  of  u  as  a  quantity  that  is  the  negative 
of  E  in  Refs  8  and  10.  We  start  by  summarizing  the  basic  algorithm  for  T  , 

T'  ,  T"  and  T"'  ,  including  such  comments  as  seem  necessary,  in  particular  in 
regard  to  the  computation  of  quantities  with  minimal  rounding  error.  The 
rounding-error  problems  referred  to  here  do  not  extend  to  the  most  serious  one 
of  all,  that  arises  for  orbitaL  paths  that  are  sufficiently  close  to  parabolic: 
the  entire  'basic  algorithm'  then  has  to  be  replaced  by  a  series-based  algorithm, 
which  i11  also  rooted  in  the  formulae  of  Lancaster  et  al  and  is  summarized  after 
the  basic  algorithm. 

The  basic  algorithm  operates  whenever  any  one  (or  more)  of  the  following 
three  conditions  applies:  (i)  m  >  0  (multi-revolution  elliptic  path);  (ii) 
x  <  0  (remote  from  parabolic  path,  for  which  x  =  1  );  (iii)  |u|  >  0.4  ,  where 


the  criterion  (0.4)  being  an  essentially  arbitrary  one.  The  appropriate  part  of 
the  basic  algorithm  also  always  operates  when  the  special  function  of  the  sub¬ 
routine  is  required,  since  there  is  never  a  need  for  series  expansion  when  the 
computation  does  not  proceed  as  far  as  T  .  The  following  formulae  apply  to 
both  cLlipses  and  hyperbolas: 


4 


y 


( B— 2  ) 


■  J-Z  o 
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z 

^(1  -  q2  +  q2x2)  , 

(B-3) 

a 

ts 

z  -  qx  , 

(B-4) 

A 

= 

z  +  qx  , 

(B-5) 

B 

■ 

qz  -  x  , 

(B-6) 

B 

= 

qz  +  x  , 

(B-7) 

f 

ay 

(B-8) 

and 

S 

= 

xz  +  qu  . 

(B-9) 

2  2 
In  equation  (B-3),  z  is  computed  using  the  input  argument  1  -  q  , 

provided  independently  of  q  ,  as  part  of  the  rounding-error  minimization  philos¬ 
ophy.  In  respect  of  the  next  four  quantities,  it  is  actually  a  and  3  ,  rather 
than  A  and  B  ,  that  are  normally  required.  Viz  a  in  (B-8)  above  and  B  in 
(B- 13)  below,  bo  the  flow  will  often  bypass  A  and  B  (required,  in  their  own 
right,  only  in  the  special  operation  of  the  subroutine).  However,  about  50%  of 
the  time  it  will  be  necessary  to  compute  a  and  B  from  A  and  B  ,  to  avoid 
rounding  error  in  subtractions,  and  sometimes  (in  the  special  operation)  the 
opposite  computations  will  be  necessary;  the  connecting  formulae  are 

aA  3  1  -  q 

and 

3B  »  (1  -  q2) (q2u  -  x2) 

It  will  also  be  necessary,  50%  of  the  time,  to  compute  g  from 

g  3  (x2  -  q2u) / (xz  -  qu) 

instead  of  from  (B-9). 

We  now  require  a  quantity,  d  ,  for  which  there  are  formulae  that  differen¬ 
tiate  between  ellipses  and  hyperbolas.  For  an  elliptic  path, 

d  “  mir  +  arg(g  ,  f)  ,  (B-12a) 


CB— 10) 
(B— 1 1 ) 
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where  the  arg  function  is  as  for  equation  (19)  of  the  main  text.  The  corre¬ 
sponding  formula  for  a  hyperbolic  path  is 

d  »  tanh  * (f/g)  ,  (B-12b) 

which  may  in  principle  be  evaluated  as  ln(f  +  g)  ,  as  assumed  by  Lancaster  et  at 
(since  f  and  g  are  equal  to,  not  just  proportional  to,  sinh  d  and  cosh  d  , 
respectively).  This  would  lead  to  serious  error  when  f  is  small,  however  (not 
noted  by  Lancaster  at  at),  so  in  this  case  we  evaluate  the  inverse  tanh  function 
by  series;  as  in  Ref  29,  it  is  more  efficient  to  evaluate  d  as 

2  tanh  1 j  £ / ( g  +  1)[  ,  rather  than  as  tanh  ' (f/g)  .  It  is  perhaps  worth  remarking 
on  what  may  appear  to  be  curious  in  the  preceding  formulae,  namely,  that  the 
elliptic  and  hyperbolic  regimes  can  be  dealt  with  by  just  two  versions  of  a  single 
formula,  equations  (B-12);  the  explanation  13  that  the  distinction  really  starts 
with  the  absolute-value  sign  in  (B-2) ,  as  a  result  of  which  (B-8)  and  (B-9)  are 

reached  wich  values  of  i  and  g  such  that  (apart  from  rounding  error) 

2  2  2  2 
g  +  f  a  1  for  ellipses  and  g  -  f  =1  for  hyperbolas. 

The  formula  for  T  can  now  be  given,  to  cover  all  orbits  that  are  not 
'too  parabolic'.  It  is 


T  -  2(d/y  +  3) /u  .  (B— 13) 

This  is  taken  directly  from  Refs  8  and  10,  where  the  formula  for  T1  Is  given 
as  well,  vis 

T'  =  (3xT  +  h q3x/z  -  4)/u  .  (B- 14) 

In  ( B— 14),  as  in  ( B— 13),  we  are  not  concerned  with  non-zero  values  of  the 
u-denominator,  since  they  only  arise  with  the  near-parabolic  rdgime  that  we  are 
not  yet  considering,  but  the  z-denominator  requires  attention.  From  (B-3)  we  see 
that  z  can  only  vanish  when  x  is  zero  and  q|  ■  1  ,  so  we  do  not  get  an 
infinity  in  (B— 14)  but  an  indeterminate  term  of  the  form  0/0  (herein  lies  a 
further  explanation  of  the  discontinuities  in  Fig  2  that  were  remarked  upon  in 
section  5.2).  The  effect  of  this  indeterminacy  is  entirely  localized,  however, 
with  no  spread  of  rounding  error  in  the  vicinity  of  the  zero  in  z  ;  thus  all 
that  is  necessary  is  that  the  computation  of  T'  be  bypassed  when  z  is  zero. 

To  complete  the  description  of  the  basic  algorithm,  it  only  remains  to 
present  the  formulae  for  T1'  and  T"'  .  These  were  not  given  by  Lancaster  et  at, 
but  are  easy  to  derive,  on  proceeding  from  (B— 14)  and  bearing  in  mind  that 
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Accuracy  demands  that  b  be  computed  with  minimal  rounding  error.  The 
n  3 

critical  computation  is  of  b^  ,  which  we  simply  set  to  1  -  q  if  q  <  0.5  . 

If  q  >  0.5  ,  however,  we  set 

bQ  ■  jq  +  1/(1  +  q)  | ( 1  -  q2)  .  (B-24) 

The  two  expressions  for  are  formally  identical,  but  computing  accuracy  (when 

needed)  is  given  by  (B-24),  because  of  the  assumed  accuracy  of  the  independent 
2 

argument  (1  -  q  )  of  TLAM13.  For  n  >  0  ,  b^  is  given  by  the  recurrence 
formula 

b  -  b  ,  +  q2n+1  (1  -  q2)  .  'B-25) 

n  n-1 

The  computation  of  is  straightforward,  since  we  may  express  it  as 

A  -  a  / (2n  +  3)  ,  (B-26) 

n  n  * 

where  an  is  given  by  the  recurrence  formula 

an  -  |(2n  -  1)/2n|an_1  ,  (B-27) 

with  =  4  , 

In  principle,  then,  T  is  computed  from  (B-21),  with  terms  continually 
added  until  there  is  no  further  change  of  value.  But  for  values  of  |u|  close 
to  the  maximum  (0.4)  for  which  the  series  would  be  used,  the  convergence  was 
found  to  be  rather  slow.  At  the  expense  of  a  little  extra  computation,  however, 

it  was  possible  to  accelerate  it  by  computing  the  series  for  Tx*"  instead  of 

2  2  .  . 

T  ,  and  then  dividing  by  x  (=  1  -  u  ,  of  minimum  value  0.84).  Writing 


« 


(B-28) 


N 

0 
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A  -  A  .  ■*  -  |(6n  +  1)/<4n2  -  1)  I  A  .  (B-30) 

ti  tv*  1  1  J  n 

Thus  the  c imputation  of  T  is  based  on  (B-28)  rather  than  (B-21). 

Two  additional  complications  arise  with  the  computation  of  T'  ,  T"  and 
T1"  .  First,  it  is  derivatives  with  respect  to  x  that  we  want,  whereas  the 
series  for  T  is  in  terms  of  u  .  Writing  T  for  dT/du  etc,  we  have 

T*  =  -2xT“  ,  (B-31) 

T"  =  -2T~  +  Ax2T*  (B-32) 

and 

T"'  =  1 2xT°  -  8x3T=  .  (B— 33) 

The  other  complication  arose  from  the  desire  to  compute  x2T  ,  T  ,  T~  ,  and  T 

efficiently  within  the  same  loop  of  code,  It  was  decided  that  negligible  accuracy 

2 

would  be  lost  tf  the  criterion  for  an  unchanging  value  of  Tx  was  used  to 
terminate  the  computation  of  all  four  quantities,  but  there  was  another  difficulty 
created  by  the  fact  that  the  initial  value  of  n  in  the  series  for  T  ,  T  and 
T~  should  be  1 ,  2  and  3  respectively;  to  overcome  this,  it  was  decided  that  the 
simplest  procedure  was  to  allow  the  initial  value  to  be  zero  in  each  case,  even 
though  initially-zero  contributions  to  T  ,  etc  would  be  computed  in  consequence. 

As  a  footnote  to  the  function  <p Cu)  ,  it  is  remarked  that  it  can  be 
expressed  in  terms  of  the  same  hypergeometric  function  as  was  involved  in  the 
quantity  X  used  in  Appendix  A.  Thus  we  have  <f> (u)  e  ^-F (3,1  ,2J;J(1  -  x))  , 
where  i ( 1  -  x)  =  !u/|l  +  /l  -  u }  when  x  >  0  ;  this  identifies  (B-18)  with 
equation  (28)  of  Ref  18.  Also,  Ref  18  indicates  a  way  of  expressing  (B-18)  with 
only  one  occurrence  of  the  function  $  ;  in  the  notation  of  the  present  Appendix, 
the  expression  for  T  is  then 

T  ”  a  |  «2d(f2)  +  4q |  . 
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FORTRAN- 7?  SUBROUTINE  TLAMB 

SUBROUTINE  TLAMB  (M,  Q,  QSQFM1,  X,  N ,  T,  DT,  D2T,  D3T) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

LOGICAL  LMI,  LI,  L2 ,  L3 

DATA  PI,  SW  /I. 14 15926 53 58 9793 DO,  0.4D0/ 

LM1  »  N.EQ.-l 
LI  ■  N.GE.  I 
L2  “  N.GE.2 
L3  »  N. EQ. 3 
QSQ  =  0*0 
XSQ  =  X*X 

U  *  (ICO  -  X) *(100  *  X) 

IF  (  .  NOT .  LM 1 )  TiiEN 

C  (NEEDED  IF  SERIES,  AND  OTHERWISE  USFFUI.  WHEN  2  - 

DT  *>  CDO 
D2T  =  ODO 
D'JT  =  ODO 
END  IF 

IF  (LM1  .OR.  M.GT.O  .OR.  X.LT.ODC  .OR.  DABS ( U) . GT . SW) 

C  DIRECT  COMPUTATION  (NOT  SERIES) 

Y  -  DSQRT ( DABS (U) ) 

Z  =  DSQRT (OSQFM1  +  QSQ*XSQ) 

QX  -  Q*X 

IF  (QX.LE.ODO)  THEN 
A  =  Z  -  QX 
B  -  Q*Z  -  X 
END  IF 

IF  (QX.LT.ODO  .AND.  LMI)  THEN 
AA  =»  QRQFMl/A 

BB  “  QSQFM1*(QGQ*U  -  XSQ)/B 
END  IF 

IF  (QX  .  EQ . ODO . AND. LMI  .OR.  QX.GT.ODO)  THEN 
AA  =  Z  +  QX 
BB  *>  Q*Z  -  X 
END  IF 

IF  IQX.GT.ODO)  THEN 
A  =  QSQFM1/AA 

B  ■  QSQFM1* (QSQ*U  -  XSQ)/BB 
END  IF 

IF  (.NOT.UU)  THEN 

IF  (QX*U.GE.ODO)  THEN 
G  =  X*Z  +  Q*U 
ELSE 

G  =  (XSQ  -  QSQ*U)/(X*Z  •  Q*U) 

END  IF 
F  ■  A* Y 

IF  (X.LE.1DO)  THEN 

T  =  !1*PI  +  DAT  AN  2  (F,  G) 

ELSE 

IF  (F.GT.SW)  THEN 
T  =  DLOG  ( F  +  G) 

ELSE 

FG1  =  F/(G  »  IDO) 

TERM  »  RD0*FG1 
FG1SQ  -  FG1*FG1 
T  =  TERM 
TWO II  =  IDO 

1  TWO I I  »  TWOI1  +  200 

TERM  =  TERM*FG1SQ 
TOLD  *  T 

T  =  T  +  TERM/TWO II 
IF  (T . NE. TOLD)  GO  TO  1 

C  (CONTINUE  LOOPING  FOR  INVERSE  TANK) 
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4  1 


END  IF 
END  IF 

T  =  2 DO* (T/Y  *  B)/U 
IF  (LI  .AND.  2.NE.0D0)  THEN 
02  »  Q/Z 
QZ2  =  Q2*Q2 
QZ  =  Q2*Q22 

DT  =  (3D0*X*T  -  4  DO*  (A  QX*QSQFM1)  /Z)  /U 

IF  (L2)  D2T  =  (3D0*T  3D0*X*DT  +  4 D0*Q2*QSQFM1 ) /l' 

IF  ( L3 }  D3T  =  (8D0*DT  *  7D0»X*D2T  -  12 DO *QZ*QZ2 '  X*CSQFM 1 1  ■ U 
END  IF 
ELSE 
DT  =  B 
D2T  -  BP 
D3T  =  AA 
END  IF 
ELSE 

C  COMPUTE  BY  SERIES 

UOI  =  IDO 
IF  (LI)  U1I  »  IDO 
IF  (L2)  U2I  -  IDO 
IF  ( L3 )  U3I  =  IDO 
TERM  =  4  DO 
TO  »  Q*QSQFM1 
1  =  0 

IF  (Q.LT.SD-1)  TOSUM  =  IDO  -  Q*QSQ 

IF  (Q.GE.SD-1)  TQSUM  »  (1D0/(1D0  *  Q)  *  Q) *QS0FM1 
TTMOLD  -  TERM/ 3  DO 
T  *  TTMOLD*TQSUM 
C  (START  OF  LOOP) 

Z  I  »  I  +  1 
P  -  I 

UOI  =  uot*u 

IF  (LI  .AND.  I.GT.l)  U1I  U1I*U 

IF  (L2  .AND.  I.GT.2)  U2I  -  U2I*U 

IF  ( L3  .AND.  I . GT . 3 )  UOI  »  U3I*U 

TERM  -  TERM* ( P  -  0.5D0)/P 

TQ  =  TQ*QSQ 

TQSUM  =  TQSUM  +  TQ 

TOLD  »  T 

TTERM  =  TERM/ ( 2D0*P  *  JDO) 

TQTERM  =■  TTERM'TQSUM 

T  »  T  -  UOI* ( ( 1 . SOO*P  -  0.25D0) ‘TQTERM/ (P*P  -  O.25D0j 
A  -  TTMOLD*TQj 

TTMOLD  -  TTFRM 
TQTERM  =  TQTERM* P 
IF  (LI)  DT  =  DT  +  TQTERM *U I I 
IF  (1,2)  D2T  =  D2T  +  TQTERM*U2 1  *  ( P  -  IDO) 

IF  ( L3 )  D3T  =  D3T  +  TQTFRM*U3 I* ( P  -  1D0)*(P  -  2DO) 

IF  (I.LT.N  .OR.  T . NE . TOLD)  GO  TO  2 
C  (END  OF  LOOP) 

IF  ( L3 )  DOT  =  BD0*X*(I.5D0*D2T  -  XSQ*D3T) 

IF  (L2)  D2T  =  2D0*(2D0*XSQ*D2T  -  DT) 

IF  (LI)  DT  =  -2DO*X*DT 
T  =  T/XSQ 
END  IF 
RETURN 
END 
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Appendix  D 

FORTRAN- 7 7  SUBROUTINE  XLAMB 


i  > 

4 


A 


SUBROUTINE  XLAMB  (M,  Q,  QSQFM1 ,  TIN,  N,  X,  /.PL; 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-2J 

PARAMETER  ( PI=3 . 14 1 59265 353979 3DO ,  TOL=3D-7,  CO=1.7DO; 

PARAMETER  (C1=0.5D0,  02=0.0300,  03=0.1500,  041=100,  C4  2 -Q . 2 4  DO. 
D8RT(X)  »  OSQRT ( DSQRT ( DSQRT ( X) ) ) 

THR2  -■  DAT  AN  2  ( QSQFM1 ,  2D0*Q)/Pt 
IF  (M.EQ.O)  THEN 

C  SINGLE-REV  STARTER  FROM  T  (AT  X  =  0]  i  BILINEAR  (USUALLf 

N  =  1 

CALL  T LAMB  (M,  Q,  Qf.QFMl,  000,  0,  TO,  OT,  D2T,  D>T) 

TOIFF  «  TIN  -  TO 
IF  (TDIFF.  I.E.ODO)  THEN 
X  =  TO*TDIFF/ ( -4D0*TIN) 

C  (-4  ts  THE  VALUE  OF  DT,  FOR  X  =  0) 

ELSE 

X  =  -TDIFF/ (TDIFF  4  4D0) 

W  =  X  4  CO*  DSQRT  ( 2D0*  ( 1  DO  -  THR2  ) ) 

IF  (W.LT.OOO) 

A  X  =  X  -  OSQRT (D3KTI-W)  )  *(X  -  OSQRT  (TDI  FF  ■' ITDI  FF  -  i  .  5  00  »  l\l  j  , 
W  =  4 DO/ ( 4 DO  *  TDIFF) 

X  *  X* ( IDO  *  X» (C1*W  -  C2*X*DSQRT(W) ) ) 

END  IF 
ELSE 

C  WITH  MULTIREVS,  FIRST  GET  T(MiN)  AS  BASIS  FOR  STARTER 

XM  »  IDO/  ( 1 . 5D0*  (M  *■  5D-1)*PI) 

IF  (THR2 . LT . 5D- 1 )  XM  =  D8RT ( 200*TWR2 ) *XM 
IF  (THR2 .GT.5D-I)  XM  =  £2D0  -  D8RT(2D0  -  2D0*THR2 , ; *XM 
C  (STARTER  FOR  TMIN) 

DO  1  I»l,  12 

CALL  TLAMB  (M,  Q,  QSQFMl,  XM,  3,  TMIN,  DT,  D2T ,  D3T) 

IF  (D2T. EQ.ODO)  GO  TO  2 
XMOLD  =  XM 

XM  *  XM  -  DT*Q2T/ (D2T-D2T  -  DT*D3T/2DO) 

XTEST  =  DABS (XMOLD/ XM  -  IDO) 

IF  (XTEST.  LE. TO L)  GO  TO  2 

1  CONTINUE 
N  =  -L 
RETURN 

C  (BREAK  OFF  5.  EXIT  IF  TMIN  NOT  LOCAIED  -  SHOULD  NEVER  HAPPE., 

C  NOW  PROCEED  FROM  T (MIN ;  TO  FULL  STARTER 

2  CONTINUE 

TDIFFM  =  TIN  -  TMIN 
IF  (TDIFFM.  I.T.  ODO)  THEN 
N  =  0 
RETURN 

C  (EXIT  IF  NO  SOLUTION  WITH  THIS  M, 

ELSE  IF  (TDIFFM, EQ.ODO!  THEN 
X  =  XM 
N  -=  1 
RETURN 

C  (EXIT  IF  UNIQUE  SOLUTION  ALREADY  FROM  X(TMIN)! 

ELSE 
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N  *  3 

IF  (D2T.EQ.0D0j  D2T  =  6D0*M*PI 

X  =  OSQRT (TDI FFM/  (  D2T/2DO  -*•  TDI FFM/  ( 1  DO  -  XM)**2Jj 
W  =  XM  +  X 

W  =  W*4D0/(4D0  +  TDI  FFM)  +  (IDO  “  W)**2 

X  *  X * ( 1 DO  -  (IDO  *  M  v  C4 1  * (THR2  -  0.5D0)J/(1D0  *  C3»M) * 
A  X*(C1*W  f  C2*X*DSQRT(W) ) )  -  XM 

D2T2  »  D2T/2DQ 
IF  (X.GE.1D0)  THEN 
N  =  I 
CO  TO  3 
END  IF 

C  (NO  FINITE  SOLUTION  WITH  X  >  XM) 

END  IF 
C 

END  tF 

C  (NOW  HAVF,  A  STARTER,  SO  PROCEED  BY  HALLEY j 

5  CONTINUE 
DO  4  1=1,3 

CALL  TLAMB  (M,  Q,  QSQFM1 ,  X,  2,  T,  DT,  D2T,  D3T) 

T  =  TIN  -  T 

IF  (DT.NE.ODO)  X  “  X  -  T*DT/(DT*DT  -  T*D2T/2Q0) 

4  CONTINUE 

IF  (N.NE.3)  RETURN 

C  (EXIT  IF  ONLY  ONE  SOLUTION,  NORMALLY  WHEN  M  ^  Oj 

N  =  2 
XPL  =  X 

C  (SECOND  MULTI-REV  STARTER) 

3  CALL  TLAMB  (M,  Q,  QSQFM1 ,  0D0,  0,  TO,  DT ,  D2T,  D3T) 

TDIFFO  =  TO  -  TMIN 
TDI FF  =  TIN  -  TO 
IF  (TDIFF.LE.O)  THEN 

X  =  XM  -  DSQRT (TDI FFM/ (D2T2  -  TDI FFM* ( D2T2/TDI FFO 
A  -  1D0/XM**2) ) ) 

ELSE 

X  =  -TDIFF/ (TDIFF  +  4DO) 

IJ  »  200 

W  -  X  +  C0*DSQRT ( 2D0* ( IDO  -  THR2 ) ) 

IF  (W.LT.ODO)  X  = 

A  X  -  DSQRT  ( D8RT  (-W)  )  *  (  X  ♦  DSQRT(TDIFF/ (TDIFF  *  !..5DO*TO))) 

W  =  400/  (4D0  i-  TDIFF) 

X  =  X* ( IDO  *  (IDO  ‘  M  -  C42* (THR2  -  0.5D0))/(1D0  -  C3*M) * 

A  X*(C1*W  -  C2»X*DSQRT(W) ) ) 

IF  (X.LE.-1D0)  THEN 
N  =  N  -  1 

C  (NO  FINITE  SOLUTION  WITH  X  <  XM) 

IF  (N. EQ. 1)  X  =  XPL 
END  IF 
END  IF 
GO  TO  5 
END 


a 
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FORTRAN~77  _SUBROUXXNE  VLAMB 


SUBROUTINE  7  LAMB  (GM,  Rl ,  R2  ,  TH,  rDELT,  N,  VR11,  7TU , 
1  VR12,  VT12,  VR21,  VT21,  VR22,  VT22) 

IMPLICIT  DOUBLE  PRECISION  <A-H,0-Zj 

PARAMETER  (PI  «  3.14159265303979)00,  TWOPt  -  2D0-PI) 

M  o  Tll/THOPI 

THR2  =  TH/ 2 DO  -  M*PI 

DR  =  Rl  -  R2 

R1R2  -  K1*H2 

R1R2TH  -  4D0*R1R2*USIN(THR2) **2 
CSQ  -  DR* *2  +  R1R2TH 
C  =  DSQRT(CSQ) 

S  *  (Rl  t  R2  f  C ) / 2 DO 
GMS  =  DSQRT (GM*S/2D0) 

QSQPM1  “  C/S 

Q  -  DSQRT ( R1R2 ) *DCOS (THR2 ) /S 
IF  (C.NE.ODO)  THEN 
RHO  =  DR/C 
SIG  »  R1R2TH/CSQ 
ELSE 

RHO  -  ODO 
SIG  =  IDO 
END  IF 

T  =  4D0*GM3*TDELT/S**2 

CALL  XLAMB  (M,  Q,  QSQFM1,  T,  N,  XI,  X2) 

PROCEED  FOR  SINGLE  SOLUTION,  OR  A  PAIR 
do  i  i=i;n 
IF  (I.EQ.l)  THEN 
X  =  XI 
ELSE 
X  =  X2 
END  IF 

CALL  TLAMB  (M,  Q,  QSQFM1,  X,  -1,  UNUSED,  Q2M1NX,  Q2PL/., 
VT2  =  GMS *ZPLQX* DSQRT (SIG) 

VR1  GMS  *  (QZMINX  -  QZPLX*RH0)/R1 
VT1  =  VT2/RI 

VR2  =  -GMS* (QZMINX  +  QZPLX*RHO) /R2 
VT2  =  VT2/R2 
IF  (I.EQ.l)  THEN 


VRI1 

- 

VRl 

VT11 

VT1 

VR12 

VR2 

VT12 

a? 

VT2 

ELSE 

VR2  1 

= 

VRl 

VT21 

ii 

VTI 

VR22 

= 

VR2 

VT22 

= 

VT2 

END  IF 

1  CONTINUE 
RETURN 
END 


UZf 
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LIST  OF  SYMBOLS 

(Principal  usage  in  the  main  text  -  some  symbols  have  different  meanings  in  the 
Appendices) 

a  semi-major  axis 

chord  of  triangle  CP^Pj 
centre  of  force 


V  P2 


T 

w,  W 


X. 

L 


*M 

+ 


eccentricity 

eccentric  anomaly 

arbitrary  non-negative  integer 

number  of  complete  revolutions  included  in  the  flight  path 

number  of  solutions  of  a  particular  Lambert  problem 

2 

semi-latus  section,  a(1  -  e  ) 

end-points  of  flight  path 

tan  i(u  -  0^)  ,  whence  also  q  »  1  -  c/s 

distances  CP1  and  CP2 

semi-perimeter  of  triangle  CP jP2 

times  such  that  At(&  t2  -  t^)  is  the  flight  time 

non-dimensionalized  At  ,  defined  by  equation  (2) 

value  of  T(x,  9)  ,  or  T(x,  q,  m)  ,  when  x  a  0 

minimum  value  of  T  for  given  0  >  2n 


M 

T' ,  eta  partial  derivatives  3T/3x  etc 
VR  radial  velocity 

transverse  velocity 

empirical  quantities  used  in  starter  weighting 

2 

fundamental  parameter  (iteration  variable)  such  that  x  *>  1  -  s/2a 
value  of  x  estimated  after  iteration  i  (x^  -  starter) 


value  of  x  associated  with  T, 


M 


second  solution  for  x  ,  when  relevant 

//.  2  2  2, 

/(I  -  q  +  q  x  ) 
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LIST  OF  SYMBOLS  (concluded) 

a  u/a 

y  /( jjs/2)  ,  introduced  at  equation  (26) 

dr,  59  absolute  values  of  residual  differences  in  r 2  and  6 

6^,  50  5r/r2  an<1 

e  smaller  of  relative  errors  in  (i)  solved-for  x  and  (ii) 

corresponding  T 

6  angle  PjCP^ 

0^  0  reduced  to  (0,  2u)  range 

0  value  of  0  for  the  Lambert-equivalent  isosceles  triangle 

\  empirical  quantity  used  in  starting  formulae 

U  strength  of  gravitational  force  centre  (C) 

p  (r1  -  r2)/c 

a  /(I  -  p2) 

x  time  from  pericentre 


Author 
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S.  Herrick 
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(doubly  degenerate  triangle) 


Fig  1  Examples  ol  Lambert-congruent  equivalent  triangles  • 
three  classes  illustrated  (same  s  ,  different  c) 
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b)  0  =  it/2 


c)  G  =  n 


d)  9  »  3n/2 


e)  0  =  2n- 


Orbital  paths  associated  with  an  evolving  isosceles  triangle 
(fixed  value  of  s  ,  with  |x|  s  0.5) 
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17.  Abstract 

During  the  past/30  years  theta  has  been  s  resurgence  of  interest  in  the  classi¬ 
cal  orbital  boundarytvalue  problem  of  Lambert,  largely  because  of  its  relevance  to 
space  rendezvous  and  interception/}  The  moat  notable  contribution  to  the  subject  was 
by  Lancaster,  Blanchard  and  Davatrfty,  in  1966,  but  more  recent  researchers  have 
failed  to  build  on  that  workr*£he  present  Report  it  aimed  at  retaadyint)  this-cegleeT*- 
Jlisy  providing  details  of  a  universal  solution  of  Lambert’s  problem, baaed  on  the 
:  approach  of  Lancaster  et  al.  <*Tn  particular, ~the~Re'port“fi’rTOBTfrk_^tarting  formulae 
for  Halley's  cubic  iteration  process,  used  for  evaluation  of  the  unknown  parameter, 
x  ,  at  the  heart  of  the  approach:  this  process  always  gives  highly  accurate  values 
of  x  after  thro*  iterations. 

A  Fortran-77  computing  procedure  for  a  general  solution  of  Lambert's  prohlem 
has  been  developed,  and  its  three  rein  subroutines  ere  listed.  Details  are  given 
of  the  testing  of  this  proce-luta. 

Much  of  who  Report  is  devoted  to  a  classification  of  the  set  of  all  Lambert 
problems,  and  to  a  discussion  of  variou  ;  groaitcic  end  physicul  aspects.  ^  , 
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